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ABSTRACT 

Recent numerical and analytical studies have shown that galaxies accrete most of their baryons 
via the cold mode, from streams with temperatures T ~ 10"* — 10^ K. At these temperatures, the 
streams should radiate primarily in the Lya line and have therefore been proposed as a model to 
power the extended, high-redshift objects known as Lya blobs, and may also be relevant for pow- 
ering a range of less luminous Lya sources. We introduce a new Lya radiative transfer code, aRT, 
and calculate the transport of the Lya emission from cold accretion in cosmological hydrodynamical 
simulations. In this paper, we describe our methodology, and address physical and numerical issues 
that are critical to making accurate predictions for the cooling luminosity, but that have been mostly 
neglected or treated simplistically so far. In particular, we highlight the importance of self-shielding 
and of properly treating sub-resolution models in numerical simulations. Most existing simulations 
do not self-consistently incorporate these effects, which can lead to order-of-magnitude errors in the 
predicted cooling luminosity. Using a combination of post-processing ionizing radiative transfer and 
re-simulation techniques, we develop an approximation to the consistent evolution of the self-shielded 
gas. We quantify the dependence of the Lya cooling luminosity on halo mass at z = 3 for the simplified 
problem of pure gas accretion embedded in the cosmic radiation background and without feedback, 
and present radiative transfer results for a particular system. While pure cooling in massive halos 
(without additional energy input from star formation and AGN) is in principle sufficient to produce 
La ^ lO'*'' — 10'*^ erg s~^ blobs, this requires including energy released in gas of density sufficient 
to form stars, but which is kept 100% gaseous in our optimistic estimates. Excluding emission from 
such dense gas yields lower luminosities by up to one to two orders of magnitude at high masses, 
making it difficult to explain the observed Lya blobs with pure cooling. Resonant scattering produces 
diffuse Lya halos, even for centrally concentrated emission, and broad double peaked line profiles. In 
particular, the emergent line widths are in general not representative of the velocity dispersion within 
galactic halos and cannot be directly used to infer host halo masses. 
Subject headings: Galaxies: formation, evolution, high-redshift - cooling flows - radiative transfer 



1. INTRODUCTION 

Steidel et al. (2000) discovered two extremely bright, 
large, and diffuse Lya-emitting 'blobs' in a narrowband 
survey of the SSA 22 proto-clustcr region at (z) — 3.09 
(for the early detection of extended Lya emission at 
z ^ 2.4, see also Francis et al. 1996, 2001; Keel et al. 
1999). These two blobs, labelled LABI and LAB2, 
have physical extent > 140 kpc, are more luminous 
(^Lya ^ 10*^ erg s~^) than typical line emitters at 
the same redshift by a factor ~ 10 — 100, and unlike 
similar halos around radio galaxies have no detectable 
radio continuum. Since their discovery, these blobs have 
become some of the most spectacular of a new class 
of sources that now counts several tens of members 
(e.g., Matsuda et al. 2004; Palunas et al. 2004; Francis 
et al. 2004; Dey et al. 2005; Nilsson et al. 2006; Saito 
et al. 2006; Smith & Jarvis 2007; Prescott et al. 2008, 
2009; Ouchi et al. 2009; Yang et al. 2009, 2010) and 
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whose nature remains unclear. Because detecting 
spatially extended Lya emission usually requires deep 
narrowband imaging, which covers thin redshift slices, 
only a small volume the Universe has been effectively 
surveyed for them to date; the blobs and their fainter 
analogues are therefore likely to be numerous and 
cosmologically significant. In fact, existing observations 
and theoretical models indicate that they may be the 
sites of massive galaxy formation and also may display 
signatures of associated active galactic nuclei (AGN) 
and/or supernova feedback. The Lya blobs thus provide 
a unique opportunity to probe the processes driving 
and regulating galaxy formation, and may be intimately 
related to phenomena including proto-clusters, mergers, 
and submillimeter galaxies. 

Although the most extreme Lya blobs have received the 
most attention, there in fact exists a wide continuum 
of spatially extended Lya sources at high redshift, for 
example the ones discovered by Saito et al. (2006) with 
line luminosities ^ 10^^ erg s~^ and those discovered 
by Ranch et al. (2008), with line luminosities as low 
as ^ 10"^^ erg s~^. Understanding the nature of these 
fainter but more numerous sources is equally important 
to develop a physical picture of galaxy formation. In 
fact, the fainter sources likely probe different (perhaps 
earlier) stages of galaxy assembly. 
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A central puzzle for the Lya blobs is that many of 
them do not appear to have a central source energetic 
enough to power their entire Lya emission (e.g., Matsuda 
et al. 2004; Nilsson et al. 2006). Even for the blobs which 
do have energetically-sufficient counterparts (e.g., Geach 
et al. 2005, 2007, 2009; Webb et al. 2009), it is unclear 
whether that energy can actually couple effectively 
to the Lya emission. Submillimeter starbursts and 
obscured AGN imply the presence of large quantities of 
dust, which acts to destroy Lya photons particularly effi- 
ciently (Neufeld 1990). It is therefore uncertain whether 
Lya photons produced by such dust-enshrouded sources 
can escape in significant amounts. Nevertheless, several 
different mechanisms been proposed to power the blobs 
and can be broadly divided into three categories: 

Embedded star formation or AGN, possibly 
obscured from direct view by dust, could photoionize 
the surrounding hydrogen nebula (Moller & Warren 
1998; Haiman & Rees 2001; Weidinger et al. 2004, 2005; 
Laursen & Sommer-Larsen 2007). 

Superv^rinds driven by starburst supernovae could 
explain the observed sizes and kinematics of the blobs, 
with the Lya emission being generated in the swept 
up material (Taniguchi & Shioya 2000; Taniguchi ct al. 
2001; Ohyama et al. 2003; Mori et al. 2004; Wilman 
et al. 2005; Geach et al. 2005). AGNs could also drive 
similar winds (e.g., Murray et al. 2005; Springel et al. 
2005; Di Matteo et al. 2005; Hopkins & Hernquist 2006). 

Cooling radiation, emitted as gas accretes onto 
forming galaxies, could produce luminous and extended 
structures. If a large fraction of the accreting gas has a 
temperature T ^ 10** — 10^ K, then most of the cooling 
radiation could be Lya (Fabian & Nulsen 1977; Katz 
& Gunn 1991; Hu 1992; Haiman et al. 2000; Fardal 
et al. 2001; Birnboim & Dekel 2003; Keres et al. 2005; 
Furlanetto et al. 2005; Dijkstra et al. 2006; Yang et al. 
2006; Dijkstra & Loeb 2009; Goerdt et al. 2010; Dayal 
et al. 2010). 

In this work, motivated by recent progress on our 
understanding of how galaxies get their gas, we focus 
on the third possibility. The methods developed could 
however be applied to the first two classes of models as 
well, and we plan to extend our calculations to model 
those processes in the future. 

In the classic sketch of galaxy formation (Rees & 
Ostriker 1977; Silk 1977; White & Rees 1978), gas falhng 
into dark matter halos is shocked and heated to the 
virial temperature. For a galaxy with a mass similar 
to that of the Milky Way, the shocked gas attains a 
temperature Tvir ^ 10^ K. In the dense inner regions of 
the halos, this gas efficiently radiates its thermal energy, 
loses its pressure support, and settles into compact 
discs where it can form stars. A wealth of recent work 
however suggests that this picture requires an important 
modification: most of the gas is never strongly shocked 
as it flows toward the central forming galaxy, but rather 
accretes in a "cold mode" , maintaining a temperature 
T < 10'"' K. Moreover, this cold accretion proceeds 
through dense filaments rather than in a spherically 



symmetric fashion. 

Although the importance of cold accretion in galaxy 
formation has only recently been demonstrated in 
high-resolution three-dimensional hydrodynamical sim- 
ulations (e.g., Fardal et al. 2001; Katz et al. 2003; Keres 
et al. 2005, 2009b; Ocvirk et al. 2008; Brooks et al. 
2009; Dekel et al. 2009), Binney (1977) had argued on 
the basis of analytic models of protogalaxy collapse 
that the amount of shock heating could be small for 
plausible physical conditions. Moreover, already in 
the first simulations of forming galaxies, most of the 
gas never heated above T ~ 3 x 10^ K (Katz & Gunn 
1991) and the importance of filamentary structures was 
recognized by Katz & White (1993) and Katz et al. 
(1994). Birnboim & Dekel (2003) carried out a stabihty 
analysis, supported by one-dimensional hydrodynamical 
simulations (see also Birnboim et al. 2007; Dekel & 
Birnboim 2006), and found that when the radiative 
cooling is efficient compared with the infall rate, the 
post-shock gas becomes unstable and cannot support 
the shock. When applied to cosmology, their results 
agree well with those of three-dimensional simulations. 

The Lya emission from cold accretion has already 
been the subject of some studies. Fardal et al. (2001) 
first evaluated the Lya cooling luminosity from hydrody- 
namical simulations and suggested that it could account 
for the Lya blobs discovered by Steidel et al. (2000). 
Haiman et al. (2000) reached a similar conclusion using 
simplified analytic arguments. Also using simulations, 
Furlanetto et al. (2005) found that the Lya cooling 
radiation from structure formation could account for 
some, but not all, of the luminosity of the Lya blobs. 
Yang et al. (2006) studied both the hydrogen and 
helium cooling radiation using simulations, but found 
the hydrogen Lya luminosity to be strongly dependent 
on the self-shielding correction applied. Recently, 
Dijkstra & Loeb (2009) developed an analytic model 
and suggested that cooling radiation from the cold mode 
could account for all the Lya blobs under reasonable 
assumptions. Using adaptive mesh refinement (AMR) 
simulations, Goerdt et al. (2010) provided supporting 
evidence for this picture. In our discussion (§5), we will 
contrast our main results with those of Goerdt et al. 
(2010), concluding that the differences with theirs most 
likely originate from the treatments of self-shielding and 
sub-resolution modeling, which are a focus of our study. 

No study focusing specifically on cooling emission 
has however combined realistic Lya radiative trans- 
fer with hydrodynamical simulations before. ^Because 
Lya photons resonantly scatter, the resultant morpholo- 
gies, spatial extents, and spectra are strongly modified 
by radiative transfer effects (e.g., Dijkstra et al. 2006). 
Since the Lya photons tend to follow paths of least resis- 
tance in space and frequency (§4), the spatial geometry 
and bulk velocity fields play critical roles in determining 
the radiation transport (for observational evidence of 

^ Other authors have included a cooling component in radia- 
tive transfer calculations of Lya-einitting galaxies (e.g., Tasitsioini 
2006a; Laursen & Sommer-Larsen 2007; Laursen et al. 2009a, a), 
but have not explicitly separated out the signatures of pure cool- 
ing or investigated the important uncertainties in detail. 
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these effects, see e.g. Kunth et al. f998; Mas-Hesse et al. 
2003) . Fully three-dimensional calculations are therefore 
necessary to make realistic predictions. Moreover, 
both the existing analytical and numerical studies 
have limitations that could induce important errors in 
quantities as basic as the integrated Lya luminosity of 
the cold streams. One such uncertainty arises from the 
exponential dependence of the Lya emissivity on the gas 
temperature for the temperatures T ^ 10^ K that are 
characteristic of cold accretion (§3). At present, most 
galaxy formation simulations do not self-consistcntly 
predict the temperature distribution within the streams. 
In fact, existing simulations usually do not follow 
the transport of the ultra-violet (UV) radiation that 
ionizes and heats dense gas. This seriously limits the 
predictive power of these calculations, since small errors 
in the temperatures can result in large errors in the 
Lya cooling emission, and potentially grossly violate 
energy conservation. As we will demonstrate, models of 
sub-resolution physics in hydrodynamical simulations 
can also introduce large errors if not properly taken 
into account. Analytical studies based on energetic 
considerations are not as sensitive to the temperature of 
the cold streams (e.g., Dijkstra & Loeb 2009), but are 
not immune of uncertainties either, since they rely on 
assumptions regarding the efhciency of Lya emission. 
Moreover, their simplified nature does not lend itself to 
detailed radiative transfer predictions. Resolving these 
issues is critical to relating the Lya emission from cold 
accretion to observations. 

Our ultimate goal is a systematic investigation of 
the Lya emission from galaxy formation that is both 
detailed in its predictions, and robust. By detailed, 
we envision predictions that can be directly compared 
with observations, and therefore require both 3D hydro- 
dynamical simulations and realistic radiative transfer. 
By robust, we mean that the predictions should be 
free of assumptions that introduce the kind of large 
uncertainties that existing studies are subject to. Due 
to the complexity of the problem, this ultimate goal is 
likely to require a long-term effort. The present paper 
is dedicated to laying down some of the foundations for 
this research program. Specifically, we present a new 
Lya radiative transfer code, named aRT, and describe 
its application to the cooling radiation in cosmological 
simulations of galaxy formation (for other applications 
of Lya radiative transfer codes to hydrodynamical 
simulations, see e.g. Cantalupo et al. 2005; Tasitsiomi 
2006a; Laursen & Sommcr-Larsen 2007; Laursen et al. 
2009a; Kollmeier et al. 2010; Zheng et al. 2010a). We 
pay particular attention to clarifying the physical and 
numerical uncertainties of these calculations, in particu- 
lar with respect to the predicted Lya luminosities, and 
illustrate the importance of radiative transfer effects. 
To do so, aside for calculating the Lya luminosities of a 
sample of halos from a cosmological volume, we focus 
our radiative transfer calculations on a particular system 
at z == 3 and explore variations in both the emission 
and radiative transfer physics. We limit ourselves to 
the most basic physical problem of accreting halos 
embedded in a cosmic ionizing background and neglect 
feedback processes. Follow up studies will build on the 
results obtained here and investigate the properties of 



the Lya emission as a function of halo mass and red- 
shift, and will extend them by incorporating additional 
physics, including feedback (Faucher-Giguere et al., in 
prep.). 

We begin by describing our hydrodynamical simu- 
lations in §2. We address the emission of Lya photons 
in §3 and explicitly demonstrate the sensitivity of 
the predicted Lya cooling luminosity on assumptions 
regarding the thermal state of the self-shielded gas. 
Using a combination of post-processing ionizing radia- 
tive transfer and re-simulation techniques, we develop 
an approximation to consistently model the evolution 
of the dense gas in the hydrodynamical simulations, 
resulting in the most robust numerical predictions to 
date. In §4, we present the results of Lya radiative 
transfer calculations for a particular system of total mass 
Mh = 2.5 X 10" M0 at z = 3 and highlight the role of 
both the bulk velocity flows and of the resonant scatters 
in shaping the emergent morphology and spectrum. 
Finally, we discuss our results and conclude in §5. The 
Appendices document the radiative transfer code aRT 
introduced in this work, the ionizing radiative transfer 
method, and relevant analytical estimates. 

Throughout, we assume a cosmology 

with (f^mj ^b7 ^Aj h, CTg, Us) = 

(0.28, 0.046, 0.72, 0.70, 0.82, 0.96), as inferred 
from the Wilkinson Microwave Anisotropy Probe 
(WMAP) five-year data in combination with baryon 
acoustic oscillations and supernovae (Komatsu et al. 
2009). While some of our hydrodynamical simulations 
were run with slightly different parameters, none of our 
conclusions are sensitive to the details of the cosmology. 
We assume hydrogen and helium mass fractions of 
X =. 0.75 and Y = 0.25 (e.g.. Buries et al. 2001), 
the coUisional ionization coefficients given in Katz 
et al. (1996), the Lya coUisional excitation coefficient 
and average number of Lya photons produced per 
recombination from Osterbrock & Ferland (2006), and 
the recombination coefficients in the appendix of Hui & 
Gnedin (1997). For convenience, some symbols used in 
this work are defined in Table 1. 

2. SIMULATIONS 

2.1. Code Details 

We compute the hydrodynamics of forming galaxies 
in a ACDM universe using a modified version of the 
GADGET cosmological simulation code (Springel 2005). 
The calculation of the gravitational force uses a combi- 
nation of the particle mesh algorithm (e.g., Hockncy & 
Eastwood 1988) for large separations and the hierarchi- 
cal tree algorithm (e.g., Barnes & Hut 1986; Hernquist 
1987) at small distances. The gas dynamics is calculated 
using a smoothed particle hydrodynamics (SPH) algo- 
rithm (e.g., Lucy 1977; Gingold & Monaghan 1977) that 
conserves both energy and entropy (Springel & Hern- 
quist 2002). The modifications with respect to the pub- 
lic version of the code include the treatment of cooling, 
the effects an uniform ultra-violet background (UVB), 
and a multiphase star formation algorithm as in Springel 
& Hernquist (2003). Star formation is implemented by 
the stochastic spawning of collisionlcss star particles by 
the gas particles. In practice, star formation in the 
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TABLE 1 

Symbols used in this work 



Symbol 



T 

Tv 

X 

Ti 

^ i,c 

CLya 

A.B, 

a. 
2:hi,3;hii 

^Hcli yHcIIi 



yHoIII 



Definition 

number density of species i 

column density of species i 

gas temperature 

optical depth at frequency v 

Lya central frequency 

Lya Doppler width 

dimensionless frequency offset (y — Vi^)l I^v-q 

photoionization rate of species i 

coUisional ionization coefficient of species % 

Lya coUisional excitation coefficient 

Lya emissivity 

case A, B recombination coefficients to species i 

fractions of hydrogen in HI, HII 

fractions of helium in Hel, HelL Helll 



multiphase model occurs above a density threshold of 
riH = 0.13 cm~^ and is calibrated to the observed Kenni- 
cutt (1998) law, although it plays only a tangential role in 
this work, which focuses the cooling emission. The ther- 
mal and ionization properties of the gas are calculated 
including all relevant processes in a plasma with primor- 
dial abundances of hydrogen and helium following Katz 
et al. (1996). 

2.2. Simulation Parameters and Halo Identification 

We use two types of simulations. To achieve high 
resolution, we 'zoom in' on individual halos within a 
larger simulation box and only follow the local gas 
dynamics at the refined resolution. This is done by first 
running a dark matter only simulation and selecting 
halos of interests. The simulation is then rerun including 
gas particles, with 8 times the original mass resolution, 
in a Lagrangian volume surrounding the halo of interest 
(e.g., Katz & White 1993). In this work, we focus on 
zoom in simulations of an individual halo (labeled Al) 
selected from a volume of side length 10 h^^ comoving 
Mpc. The Al halo has a total mass 2.5 x 10" M© at 
z = 3. As it is also important to understand the trends 
and variance between different halos, we simulated an 
entire cosmological volume consisting of a cubical box 
with a side length of 40 h~^ comoving Mpc. While the 
resolution in this volume is more limited by computa- 
tional constraints, it provides us with a large number of 
halos of different masses. Table 2 lists the simulations 
used in this work and their parameters. The minimum 
gas smoothing length is set to 0.1 of the gravitational 
softening in all the simulations. Given the importance of 
self-shielding (§3), we rerun our simulations with exactly 
the same parameters, but with the UVB artificially 
turned off in regions exceeding a certain density (suffix 
_ssUV), to be discussed in §3.2. All the simulations are 
also rerun with the star formation model turned off. 
The simulations without star formation are identified by 
the additional suffix _noSF. The Al zoom in simulations 
assume a variant of the Haardt & Madau (1996) model 
of the UVB, while our cosmological simulations use the 
more recent model of Faucher-Giguere et al. (2009). 

A fricnds-of- friends (FoF) algorithm (e.g., Davis 
et al. 1985) with linking length set to 6 = 0.2 in units 
of the mean intcrparticle separation is used to identify 
the dark matter halos in the simulations. The total 



mass of the particles within each FoF group, MpoP, 
corresponds approximately to the mass in a sphere of 
mean interior density 180 times the background matter 
density, Migob (White 2002). We therefore define 
the virial radius of a halo as the radius of a sphere 
containing Misob, r^ir = ?'i806 « [MFoF/2407rp„(2;)]i/3, 
where /9„(z) = Pcrit^m(l + z)^ and pcrit is the critical 
density at 2; = 0. In what follows, we use M as a 
shorthand for MpoF ~ Afisofc- The center of each halo 
is determined as the point deepest in the gravitational 
potential. Since we run the FoF algorithm on the dark 
matter particles only, we multiply the returned masses 
by r2„i/(i7m — Qh) to account for the baryons. 

As we want to isolate the cold accretion cooling 
radiation, the simulations studied in this work do no 
include galactic winds or AGN feedback. It is of course 
likely that the results would be somewhat modified 
if these processes were included. We plan to inves- 
tigate the effects of feedback and their observational 
manifestations in future work. 

2.3. Ionization and Thermal Structure 

Our basic SPH simulations, like most cosmological 
simulations to date, assume a UVB that is spatially 
uniform throughout the simulation volume and calculate 
the local ionization state of the gas assuming photoion- 
ization equilibrium with this background. This approach 
misses the effects of self-shielding: where the optical 
depth to ionizing photons from the background is of 
order unity or more, the gas would in reality be exposed 
to an attenuated ionizing field. This has consequences 
for both the ionization state and the thermal properties 
of the gas. Indirectly, the gas dynamics is also affected 
by the modification of pressure forces. 

The omission of self-shielding implies that (in the 
absence of local sources) dense gas in the simulations 
sees an ionizing flux stronger than in reality. As a result, 
the gas tends to be overionized. This is important for 
the Lya radiative transfer problem for two reasons. 
First, the Lya emission mechanisms which seed the 
Lya photons depend not on the total gas density, but on 
the number densities of ions (§3.1). Second, the trans- 
port of Lya photons depends on the neutral hydrogen 
distribution, as only this ion provides scattering opacity 
(§4). 
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TABLE 2 
Hydro DYNAMICAL Simulations 



Name 

zl0nl28_Al'* 
zl0nl28_Al_noSF 
zl0nl28_Al_ssUV 
zl0nl28_Al_ssUV_noSF 



L (/i-i Mpc)=» 

10 
10 
10 
10 



2 X 256^ 
2 X 256^ 
2 X 256^ 
2 X 256^ 



0.8 
0.8 
0.8 
0.8 



Description 

8x mass refinement zoom 

no star formation 

UVB off at nn > 0.01 cm"^ 

UVB off at nil > 0.01 cm^^, no SF 



gdm40n512'= 


40 


2 X 512-=' 


1.6 


full box 


gdm40n512_noSF 


40 


2 X 512^ 


1.6 


no star formation 


gdm40n512_ssUV 


40 


2 X 512^ 


1.6 


UVB off at riH > 0.01 cm^^ 


gdm40n512_ssUV_noSF 


40 


2 X 5123 


1.6 


UVB off at nu > 0.01 cm-^, no SF 



^Comoving box side length. 

''Effective total number of dark mattcr+gas particles in the box, after zoom refinement. 

•^Comoving Plummer equivalent gravitational softenini 

''Dark matter, gas, and stellar particle masses are 4 X 10^ h 
respectively. 

''Dark matter, gas, and stellar particle masses are 3 X 10^ h^' Mq, 6 X 10^ h~^ Mq, and 3 X 10^ h~^ M, 
respectively. 



length, after zoom refinement. 

' Mq, 8 X 10^ /i-i Mq, and 4 X 10^ h'^ Mq, 



Self-shielding also has important effects on the thermal 
evolution of the gas. As Fardal et al. (2001) pointed 
out, neglecting self-shielding in a simulation with a 
prescribed uniform UVB introduces an artificially high 
rate of photoheating in dense regions and thus results in 
overestimated temperatures in these regions. Moreover, 
the presence of a penetrating ionizing background 
suppresses the cooling function (e.g., Efstathiou 1992; 
Katz et al. 1996; Weinberg et al. 1997; Keres et al. 2005; 
Wiersma et al. 2009) and these effects could be amplified 
by the lack of a dynamical response of the gas to cooling, 
which would tend to make it denser and hence to cool 
even more rapidly. Although these effects are critical to 
accurately predicting the Lya cooling luminosity (§3), 
previous Lya studies have either neglected them, made 
simplifying but not necessarily correct assumptions 
(Fardal et al. 2001; Goerdt et al. 2010), or have explored 
a range of prescriptions (e.g., Furlanetto et al. 2005; 
Yang et al. 2006). 

We improve significantly over previous work by 
performing ionizing radiative transfer in post-processing 
to identify the self-shielded gas, rather than rely- 
ing on simplified criteria (§3.2). Since knowing the 
distribution of neutral hydrogen is a fundamental 
component of Lya radiative transfer, post-processing 
ionizing radiative transfer has previously been used 
by other groups for related problems (e.g., Cantalupo 
et al. 2005; Laursen et al. 2009a; Kollmeier et al. 
2010; Zheng et al. 2010a), but never before in focused 
studies of cooling emission. As we will show, the 
predicted Lya cooling luminosity is very sensitive to 
the state of the self-shielded gas. In order to obtain 
more robust predictions, we rerun our hydrodynamical 
simulations with the ionizing background turned off in 
regions above a certain density threshold (informed by 
our ionizing radiative transfer calculations) as an ap- 
proximation to the consistent treatment of self-shielding. 

For the moment, we pause to discuss the star- forming 
gas, which also requires special treatment. 



2.4. The Multiphase ISM 



The star-forming gas particles in the multiphase model 
carry effective ionic densities and temperatures that are 
mass-weighted averages of the hot and cold components 
(Springel & Hernquist 2003). While most of the mass in 
this model is in the T — 1,000 K cold component,* the 
hot component is in general much hotter (T = 10^ — 10* 
K). This results in high effective temperatures with 
simultaneously large neutral fractions, and therefore 
in high coUisional excitation rates (§3.1.2). As we 
will show, using the effective multiphase temperatures 
and ionic densities for the star-forming particles yields 
artificially high cooling luminosities, owing to the 
non-linearity of the emissivity function with respect 
to density and temperature. Moreover, in the multi- 
phase model, supernovae are responsible for pressurizing 
the ISM and are therefore an additional source of energy. 

To obtain realistic results uncontaminated by feed- 
back energy, it is necessary to exclude the star-forming 
particles from the cooling luminosity calculations. We 
explore two ways of doing this. First, we use the 
simulations with star formation, but ignore all the 
multiphase particles in the luminosity calculation. For 
this case, we assume that the star-forming particles are 
effectively optically thin for the purpose of transporting 
the Lya photons. In reality, dust may destroy a portion 
of those photons, but we do not model this effect for 
this extreme prescription. This case approximates 
a multiphase medium in which cold neutral clumps 
embedded in a hot medium are either so compact that 
their covering factor is negligible, or in which they 
simply reflect the Lya photons and contribute only a 
small effective Lya optical depth (e.g., Neufeld 1991; 
Hansen & Oh 2006). A potential worry with this 
approach, however, is that it misses the cooling that 
occurs in particles with density above the star formation 
threshold. In our second approach, we make sure to 
capture all the cooling by using identical simulations 
but with the star formation model turned off, in which 

* The terminology with respect to temperature is somewhat dis- 
crepant in the contexts of galaxy formation and of the interstellar 
medium (ISM). While the T = W^ K gas is termed cold in galaxy 
formation and throughout most of this paper, it is usually qualified 
as warm in the context of the ISM to distinguish it from the much 
cooler, star-forming molecular gas. 
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Fig. 1. — Lya einissivity per n^. Left: Gas in ionization equilibrium with a hydrogen photoionization rate Fhi = 10^^^ s^^. Right: 
Gas in pure coUisional ionization equilibrium. The solid curves show the Lyo emissivity from collisional excitation and the dashed curves 
show the emissivity from recombination. For the photoionization equilibrium case, the collisional excitation curve depends on njj; we show 
log (nn/cm"^) = —6, — 4, and — 2 with curves of increasing thickness. The emissivity from collisional excitation scales with rim °^ ^ui 
in the photoionization equilibrium case. The unavoidable effects of collisional ionization are included in the photoionization equilibrium 
case. These curves assume that all the helium is in the form of Helll. 

case we can simply sum over all the particles. The gas is 
allowed to become arbitrarily dense (as permitted by the 
resolution) and the absence of ionizing and mechnical 
feedback from embedded stars makes the density peaks 
very optically thick to Lya photons. By comparison 
with the case of optically thin star-forming gas, this 
prescription therefore allows us to also quantify the 
effects of the opacity provided by star-forming particles 
on the transport problem. 

We intentionally do not model the Lya photons 
produced by stars in this work in order to separate 
out the properties of pure cooling emission; we briefly 
discuss their importance in the discussion (§5) and in 
Appendix A. 

3. lyq emission 

3.1. Emission Processes 

3.1.1. Recombination 

Ionizing radiation (either from the cosmic background, 
local star formation, or an AGN) can photoionize gas 
that recombines and produces Lya photons. Collisions in 
gas of sufficient density and temperature can also ionize 
hydrogen and be followed by the recmission of Lya pho- 
tons via recombination. We group these two processes in 
'recombination emission.' Recombination emission pro- 
duces Lya photons at a rate (in units of ph s~^ cm~^) 



ph,rcc 



= .fa,-[ccaHi(.T)nmine 



(1) 

where a^j(r) oc T^o'^ is the hydrogen case B recom- 
bination coefficient and /a,roc is the average number of 
Lya photons produced per case B recombination. For gas 
at r = 10'' K that is optically thick to Lyman series tran- 
sitions (so that higher-order Lyman series recombination 
photons can ultimately be degraded into a Lya photon) , 
/a,rec = 0.68 (Ostcrbrock & Ferland 2006). This frac- 
tion is only weakly dependent on temperature and so we 
assume this constant value throughout. 



3.1.2. Collisional Excitation 

Collisions can also excite the Lya line without ionizing 
hydrogen. The Lya emissivity from this process is given 

by 

ePh.-ii = CLya(r)nHine, (2) 

where Ci^yaiT) ex T~^/^ exp {—hva/kT) is the Lya colli- 
sional excitation coefficient in units of ph cm'^ s~^. Note 
that the Lya recombination and collisional excitation 
emissivities scale differently with temperature. More- 
over, whereas the recombination term is proportional to 
the HII number density, the collisional excitation term 
is proportional to the HI number density. The relative 
importance of the two processes will therefore depend on 
the local temperature and ionization state of the gas. 

3.1.3. Limiting Equilibrium Cases 

We assume ionization equilibrium, which is generally 
vahd since the time scale tcq = [Fhi + (Fhlc(7") + 
aj|j(T))ne]~^ to reach equilibrium is small compared to 
the dynamical time scale in both optically thin and self- 
shielded gas. The statistical equilibrium equation for hy- 
drogen is 

Fhi"hi + rHi,c(T)nenHi = aHi(T)nenHii, (3) 

where Fhi is the photoionization rate, FHi,c(r) is the 
collisional ionization coefficient, and a^j(r) is the case 
A recombination coefficient. 

Physical intuition can be gained by considering the 
two limiting cases of photoionization equilibrium 
and of pure collisional ionization equilibrium (GIF; 
Fhi ^ rHi,c(r)ne). As we will show, the two cases are 
directly relevant to our problem: While most of the cos- 
mic volume is well approximated by the photoionization 
equilibrium regime, the dense cold gas (including the 
cold streams of interest) can self-shield from the external 
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TABLE 3 
Prescriptions for Calculating the Lyo Cooling Luminosity 



# Description 

1 Standard hydro (uniform UVB and multiphase SF model), sum all particles 

2 Standard hydro, with post-processing ionizing RT, no Lya from self-shielded gas 

3 Standard hydro, with post-processing ionizing RT, CIE with Tcie = 10, 000 K 

4 Standard hydro, with post-processing ionizing RT, CIE with Tcie = 15, 000 K 

5 Hydro with uniform UVB but no SF, sum all particles 



Notes 

Overestimates Lyo luminosity 

due to UVB and multiphase model 
Satisfies energetic bound 
Ri same Lya luminosity as 2 
Ri lOx more Lya luminous than 3 
Overestimates Ly« luminosity 
due to UVB only 



Self-shielding approx. for the UVB, with SF, sum all particles 

Self-shielding approx. for the UVB, no SF, sum all particles 

Self-shielding approx. for the UVB, no SF, only sum nn < 0.13 cm~^ 

Self-shielding approx. for the UVB, SF excluded in post-processing, sum all particles 



Overestimates Lya luminosity 

due to multiphase model only 
Consistent Lyo luminosity estimate 
Consistent Lyo luminosity estimate 
Consistent Lyo luminosity estimate 
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Fig. 2. — Lya luminosity within the virial radius as a function of halo mass at z = 3 calculated from our cosmological volume simulations, 
for different physical assumptions. Each point corresponds to a randomly selected halo. The different cases are defined in Table 3, where 
corresponding remarks are given. Le/i.- Standard hydrodynamical simulation with uniform UVB and a multiphase star formation model, in 
some cases post-processed with the ionizing radiative transfer scheme to identify the self-shielded gas, corresponding to prescriptions 1, 2, 3, 
and 4. The green -|-s show the result of a simulation with a uniform UVB but no star formation, where the cooling luminosity is integrated 
over all the particles (prescription 5). Right: Simulations with the ionizing background turned off in regions where njj > 0.01 cm~^ to 
approximate self-shielding, corresponding to prescriptions 6, 7, 8, and 9. For the last three (most consistent) cases, the filled symbols show 
the results for the Al halo in our zoom in simulation with 8x better mass resolution. Caution should be exercised when interpreting the 
quantitative details of the lowest-mass halos shown, as the hydrodynamics may not be fully converged (§3.2.2). The dashed lines show an 
analytic estimate for the maximum average cooling luminosity available from the release of gravitational potential energy (Appendix A); 
the dotted lines show the more sophisticated analytic model of Dijkstra & Loeb (2009) for their fiducial parameter /grav = 0.3. 



ionizing radiation and is then more accurately described 
by the the pm^c CIE case. Figure 1 shows £^/n'^ from 
both recombination and coUisional excitation. The left 
panel shows the case of gas in ionization equilibrium with 
a photoionization rate Fhi = 10^^^ s^^ (approximately 
the magnitude of the cosmic ionizing background at 
z « 2 — 4, e.g. Faucher-Giguere et al. 2008a,b), and the 
right panel shows the case of gas in coUisional ionization 
equilibrium. The most important point to note is that 
the coUisional excitation contribution is exponentially 
sensitive to temperature, at the temperatures T ~ 10** 
K characteristic of cold accretion streams, in both cases. 
As coUisional excitation is the dominant Lya cooling 
mechanism, it is critical to accurately capture the 
thermal state of the emitting gas. The curves shown in 
Figure 1 assume that all the helium is in the form of 
Helll for simplicity. The emissivity from both processes 
is only weakly sensitive to the helium ionized fractions, 



except for temperatures T < 10'* K for the CIE case and 
at very high densities in the photoionization case, when 
most of the hydrogen is neutral and helium dominates 
the free electrons. Since little Lya emission originates 
from these regimes, our results are broadly insensitive to 
the helium ionization state, although we do solve for the 
correct helium ionization fraction in self-shielded regions 
in our simulations with the on-the-fly self-shielding 
approximation (§3.2.2). 

In our simulations, the Lya luminosity is evaluated di- 
rectly from the SPH particles. Specifically, each particle 
is assigned a Lya luminosity L^^^ = Vp{eP^'''"' + eP'^'™*), 
where the emissivity terms are evaluated using its den- 
sity, ionization state, and temperature, and Vp = Mp/pp 
is its volume, defined as the ratio of its mass to its 
density. This approach is desirable as it accurately takes 
into account the clumping of the gas on small scales, rel- 
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Fig. 3. — Hydrogen neutral fraction vs. total hydrogen number density in a cube of side length 1 comoving Mpc/h centered on the Al 
system at ^ = 3. Left: Values for a standard simulation with a uniform UVB. Right: Same quantity after post-processing with the ionizing 
radiative transfer method. The ionizing radiative transfer shows that the main effect of self-shielding is to create a vertical "plume" of 
neutral gas above a density njj ^ 0.01 cm"''. Figures 4 and 5 illustrate the effects of self-shielding on the dynamics and thermal state of 
the gas, when it is approximated during the course of the hydrodynamical simulation. The 2D histograms are in arbitrary (but matching) 
logarithmic units and weighted by n^ to emphasize the regions where the two-body emission processes are most efficient. Gas from 
multiphase, star-forming particles (njj > 0.13 cm"'') is excluded. 



evant to calculate the emission from the density-squared 
processes, and because it avoids artificial mixing that 
could occur if hot and cold phases were averaged in a 
gridding procedure. Such artificial mixing could boost 
the predicted Lya luminosity by a large factor owing to 
the non-linearity of the emissivity function. 

3.2. Self-Shielding 

Since the emission processes scale with density squared 
(eqs 1-2), the emissivity peaks in the densest regions, 
which are the most likely to self-shield. Because our 
hydrodynamical simulations lack proper ionizing radia- 
tive transfer, they do not correctly capture self-shielding 
(§2.3). As explained in §2.4, naively integrating over 
multiphase SPH particles could also induce large errors 
in the predicted cooling luminosity. As we will show, it 
is necessary to both exclude multiphase particles from 
the calculation and to model self- shielding to accurately 
predict the cooling luminosity. The rest of this section 
is dedicated to demonstrating the importance of each 
potential source of error and to developing a consistent 
approximation to the Lya cooling luminosity. 

We follow the following steps: 

1. Naively calculate the Lya luminosity from a simu- 
lation with standard UV background and star for- 
mation treatments. 

2. Using post-processing ionizing radiative transfer, 
identify the self-shielded gas in step 1. 

3. Using the post-processed output, illustrate how the 
predicted Lya luminosity depends on the assumed 
thermal state of the self-shielded gas. 

4. Rerun a hydrodynamical simulation with the same 
initial conditions, but with the ionizing background 



turned off in self-shielded regions on the fiy as an 
approximation to the self-consistent effects of self- 
shielding. 

5. By rerunning identical simulations with star forma- 
tion turned off, separate the effects of incorrectly 
including multiphase particles from those of ignor- 
ing self-shielding. 

3.2.1. Post-Processing Self- Shielding 

A technical description of our ionizing radiative 
transfer code is provided in Appendix B. Briefiy, the 
hydrogen photoionization rate of the cosmic background, 
Fjjj^, is specified and taken as the boundary condition 
at the faces of the cubical radiative transfer volume. 
For the radiative transfer calculations, the simulation 
outputs are interpolated onto a Cartesian grid taking 
into account the smoothing kernels, with A^p grid points 
along each dimension. We employ the fiducial choice 
A'p = 256 and a radiative transfer volume of (1 comoving 
Mpc/h)'^, centered around each halo considered, which 
convergence tests suggest is sufficient (§4.3). Rays 
normal to each of the six faces are then sent inward and 
the optical depth to ionizing photons is calculated along 
each ray. Given the attenuated photoionization rate at 
each point, the ionization equilibrium is updated taking 
into account photoionization, collisional ionization, 
and recombination. The procedure is iterated until 
the ionized fraction has converged in all the cells. In 
solving for the equilibrium ionization balance, the gas 
temperatures used are those provided by the hydrody- 
namical simulation. These should be accurate in the 
optically thin regions and therefore our scheme should 
accurately capture the onset of self-shielding. In this 
post-processing treatment, the temperature structure 
will however be inaccurate in the self-shielded regions. 
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Fig. 4. — Hydrodynamical properties of the Al system at z = 3 as a function of the self-shielding treatment. The left column shows 
the total projected gas mass, the central column shows the neutral hydrogen column density, and the right column shows the projected 
gas temperature. In all cases, the projected depth is 1 comoving Mpc/h. The projected temperature is weighted by density squared to 
emphasize the dense filaments and the virial radius of the halo is indicated by the dashed circles. Top: Standard simulation output, with 
a uniform ionizing background. Middle: Same, but post-processed with the ionizing radiative transfer scheme to identify the self-shielded 
regions (§3.2.1). Bottom: Same initial conditions but with the on-the-fly self-shielding approximation, in which the ionizing background is 
turned off in regions with rajj > 0.01 cm~^ as the simulation proceeds to capture the effects on the thermal and dynamical evolution of the 
gas. Figure 5 illustrates the effects of self-shielding on the temperature structure of the gas more explicitly. The simulations shown include 
star formation. 



since the modifications of the heating and coohng func- 
tions are not properly modeled in the hydrodynamical 
calculation. As outlined above, we address this in two 
ways: first, we explore a range of prescriptions for 
the self-shielded gas, illustrating the sensitivity of the 
predictions to these prescriptions; we then subsequently 
improve the accuracy of our calculations by approxi- 
mating the self-consistent thermal evolution of the gas 
with simulations in which the ionizing background is 
switched off in dense regions. The self-shielded cells 
are defined as those that see an attenuated ionizing 
background, (e"^'> < (e"'^')crit with (e"'^')crit = 0.1, 
where (e""^') is the angle-averaged attenuation factor in 
the cell after post-processing. Since the optical depth 
rapidly increases within a self-shielded region, the results 
are weakly sensitive to the choice of the self-shielding 



threshold. 

To investigate how the Lya luminosity depends on 
assumptions regarding the self-shielded gas (steps 1, 2, 
and 3 above) , we start with hydrodynamical simulations 
with standard treatments of the UV background and 
of star formation. Halos covering a broad range of 
masses are selected from our cosmological volume 
simulation gdm40n512 at z = 3 and the luminosity of 
each is defined as the sum of the luminosities of the 
gas particles contained within its virial radius. The 
same prescriptions are also applied to the Al halo 
at z = 3, which will be used for the Lya radiative 
transfer calculations (§4). We explore three cases that 
are illustrative of the range of possibilities (for similar 
prescriptions, see Furlanetto et al. 2005): 
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Fig. 5. — Gas temperature vs. total hydrogen number density in a cube of side length 1 comoving Mpc/h centered on the Al system at 
2 = 3. Left: Values for a standard simulation with a uniform UVB (top panel of Fig. 4). These also apply for the simulation post-processed 
with the ionizing radiative transfer scheme, since it does not update the temperatures (middle panel of Fig. 4). Right: Same quantity, 
but for the simulation with the ionizing background turned off in regions with nn > 0.01 cm~^ during the course of the hydrodynamical 
calculation as an approximation to the effects of self-shielding (bottom panel of Fig. 4). The self-shielded gas is generally cooler (with 
T < lO** K) when its evolution is consistently modeled as a result of the suppression of artificial photoheating and the enhancement of its 
cooling function. The 2D histograms are in arbitrary (but matching) logarithmic units and weighted by n^ to emphasize the regions where 
the two-body emission processes are most efficient. Gas from multiphase, star-forming particles (njj > 0.13 cm~^) is excluded. 

gas cools very effectively, so that higher temperatures 
are not expected. Furthermore, the proximity of 15,000 
K to the peak of the cooling curve implies that this case 
is already quite optimistic. 



Naively sum all the particles. In this simplis- 
tic prescription, we simply sum all the SPH particles 
within the virial radius. This includes dense particles 
that would in reality self-shield but that are artificially 
illuminated by a uniform ionizing background, and 
star-forming particles that carry effective multiphase 
values for their temperature and ionization state, which 
will introduce luminosity errors as described in §2.4. 

No Lya emission from self-shielded gas. Self- 
shielded gas may become neutral and cool substantially. 
If the gas is not photoionized and cools below T w lO"' 
K, its Lya emissivity is severely suppressed (Fig. 1). To 
illustrate how this case might differ from the naive cal- 
culation above, we model it by the extreme assumption 
that self-shielded gas does not produce Lya photons at 
all. To transport the Lya photons, we assume that this 
gas has a temperature T = 10, 000 K. 

Collisional ionization equilibrium. Somewhat 
intermediate between the two above cases, self-shielded 
gas could settle to collisional ionization equilibrium 
(CIE) with a temperature Tcie ^ 10, 000 K if gravita- 
tional heating is sufficiently efficient. In this case, the 
neutral hydrogen fraction is given by 

""^ " l + rHi,c(TciE)/a^i(TciE) ^^^ 

and the Lya emissivity ea/n^ is a well defined function 
of temperature, as shown in the right panel of Figure 
1. We explore how the luminosity depends on the 
prescribed CIE temperature, for Tcie = 10, 000 K 
and 15,000 K. As we will see in our simulations that 
approximate the effects of self-shielding, the self-shielded 



For the latter two cases, in which we assume ei- 
ther no emission from self-shielded gas or CIE, we 
exclude emission from the star-forming particles that 
might fall outside of the self-shielded regions, in order 
to avoid potential confusion with artificially high lumi- 
nosities from multiphase particles. Table 3 summarizes 
the different prescriptions explored for calculating the 
Lya cooling luminosity; the above cases are labeled 1—4. 

The left panel of Figure 2 shows how the Lya cooling 
luminosity varies with halo mass at the fiducial redshift 
2; = 3 for these different prescriptions. The dashed 
curve shows an analytic estimate of the average max- 
imum cooling luminosity achievable from the release 
of gravitational potential energy as a function of halo 
mass (see Appendix A; also Goerdt et al. 2010 for 
similar ideas). Briefly, the ACDM cosmology predicts 
the average mass accretion rate onto dark matter halos 
as a function of mass and redshift (e.g., Neistein & 
Dekel 2008; McBride et al. 2009; Fakhouri et al. 2010), 
as well as the shape of the dark matter halo potential 
wells (e.g., Hernquist 1990; Navarro et al. 1997). The 
product of the halo potential well depth with the gas 
mass accretion rate provides an estimate of the rate 
at which gravitational potential energy that can be 
radiated is "injected" into the halo. Figure 2 assumes 
an optimistically high efficiency factor /off = 1 (eq. 
A3). As we discuss in Appendix 2, this gravitational 
power estimate is not strictly an upper bound for the 
total amount of cooling achievable even in the case of 
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pure accretion, without feedback from stars and AGN, 
since the accretion streams are embedded in the cosmic 
ionizing background, which can transfer further energy 
to them. This contribution is however included in our 
simulations and circumstantial evidence suggests that it 
does not dominate. In fact, in our most realistic simu- 
lations with a consistent self-shielding approximation - 
shown in the right panel of Figure 2 and to be discussed 
below - the gravitational power upper bound is never 
systematically violated, and a simulation of our Al 
halo with the ionizing background completely turned 
off yields a nearly equal cooling luminosity (within 
30%) at 2; = 3 as our prescription 9 when excluding the 
multiphase star forming regions in the same manner. 

The more sophisticated analytic model of Dijkstra 
& Loeb (2009), shown by the dashed lines for their 
fiducial efficiency parameter /g^av = 0.3, includes a 
factor fcoid{M) accounting for the decreasing fraction of 
cold gas in massive halos (e.g., Keres et al. 2005, 2009b). 
This model therefore predicts lower cooling luminosities 
than the above upper bound, with a shallower mass 
dependence at large masses that is in better agreement 
with the simulation data points. It is important to 
note here that the cooling luminosities shown in Figure 
2 are "theoretical" or "intrinsic" , meaning that they 
include all the photons emitted within the virial radii 
of the halos. These luminosities will in general be 
higher than the observationally inferred luminosities, 
which only include the emission above a certain surface 
brightness threshold determined by the observation. 
Furthermore, a certain fraction of the emitted photons 
are in practice absorbed by the intervening intcrgalactic 
medium (IGM). The Dijkstra & Loeb (2009) data 
points plotted here (which were computed for a 50% 
IGM transmission factor by these authors) have been 
multiplied by a factor of 2 for a fair comparison with our 
simulation data points (which assume 100% transmis- 
sion) . In future work, we will quantify how the predicted 
theoretical luminosities translate into observational ones. 

There are two main points to take away from the 
left panel of Figure 2. First, the predicted Lya coohng 
luminosity is extremely sensitive to the treatment of 
the dense gas, with the range exceeding three orders of 
magnitude at Mh ^ 10^" M0. Second, some prescrip- 
tions actually lead to unphysical results, as comparison 
with the analytic upper bound indicates that they emit 
more Lya power than is available from the release 
of gravitational energy by orders of magnitude (in 
Appendix A, we show that photoionization from the 
cosmic background cannot physically account for such 
large luminosities either). This is the case, in particular, 
for a standard simulation with optically thin ionizing 
balance and a multiphase star formation model in 
which all the gas particles are naively summed over (red 
x's) . Unsurprisingly, the cases of no Lya emission from 
self-shielded gas (magenta squares) and of CIE with 
TciE = 10, 000 K (blue circles) yield nearly identical 
luminosities, since the cooling curve is already strongly 
suppressed at this temperature. The case of CIE with 
2ciE — 15, 000 (cyan diamonds) yields more optimistic 
Lya luminosities, although these push the limit of the 
power that can be provided by the release of gravita- 



tional potential energy alone (the dashed line in Fig. 
2) at low masses. Prescription number 5 (green -l-'s), 
which uses a simulation in which the multiphase star 
formation model was turned off, but with a uniform 
ionizing background penetrating deep into the dense 
gas, illustrates that artificial photoheating alone can 
boost the cooling luminosities by orders of magnitude. 

The critical question is therefore: What is the cor- 
rect Lya cooling luminosity of the cold streams? To 
address this, we develop a simple approximation to 
the self-consistent evolution of the gas properties with 
self-shielding. When self-shielding is properly modeled, 
it will also be possible to show that naively integrating 
over multiphase particles alone can also artificially boost 
the cooling luminosity by a large factor. 

3.2.2. On-the-Fly Self- Shielding 

In Figure 3 we plot the hydrogen neutral fraction 
a;Hi ^ JT-Hi/^H as a function of total proper hydrogen 
number density nn for the gas around the Al halo at 
z — i. The panel on the left shows this distribution for 
the standard simulation with a uniform ionizing back- 
ground. The panel on the right shows exactly the same 
quantity after the post-processing ionizing radiative 
transfer. The effect of self-shielding is clear. Roughly, 
it generates a vertical "plume" above nn ~ 0.01 cm~^, 
indicating the fact that the gas becomes mostly neutral 
above this density. This motivates our approximation 
to the self-consistent evolution of self-shielded gas in 
the hydrodynamical simulations. Namely, we rerun 
simulations with exactly the same initial conditions 
and other physical parameters, but set the ionizing 
background to zero in regions where the density exceeds 
the fiducial threshold nn = 0.01 cm""^ (for an alter- 
native scheme in which the UVB is turned off where 
the gas is optically thick to ionizing photons on a scale 
^ 0.1 — 1 kpc, see Sommer-Larsen 2006; Laursen & 
Sommer-Larsen 2007; Laursen et al. 2009b). By turning 
the ionizing background off in dense regions on the fly, 
their thermal and dynamical properties are consistently 
evolved with the modified cooling and heating functions, 
and the corresponding dynamical response. The "CDB" 
simulation analyzed by Goerdt et al. (2010) employed 
an analogous scheme, but with a density threshold 
10 X higher, nn = 0.1 cm~^; in §5, we argue that this 
difference likely explains much of the discrepancy with 
our results. Some uncertainty is introduced by our 
choice of a fixed density threshold for self-shielding, 
and in the future it would be useful to improve the 
methodology by performing proper ionizing radiative 
transfer on the fly, which our codes do not allow us to 
do at present. There are however reasons to believe that 
this choice is a good one, which we outline next. 

Figure 4 summarizes the hydrodynamical proper- 
ties (total gas distribution, neutral gas distribution, 
and temperature structure) for the Al system at z = 3 
for the different treatments of self-shielding: standard 
uniform ionizing background, post-processing ionizing 
radiative transfer, and the on-the-fly self-shielding 
approximation. When the ionizing radiative transfer is 
taken into account, the neutral hydrogen column density 
of the cold streams can be greatly enhanced, especially 
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in the higher density regions close to the central and 
satellite galaxies, indicating the fact that they self-shield 
(some of the cold gas at larger radii however remains 
optically thin). The ionization structure obtained with 
the on-the-fly self-shielding approximation is further- 
more remarkably similar to the one obtained with the 
post-processing ray tracing scheme, supporting the 
validity of using our simple density criterion during the 
course of the hydrodynamical simulation. The simula- 
tion with on-the-fly self-shielding is the most accurate 
as it consistently captures the thermal evolution and 
dynamical response of the self-shielded gas. Figure 5 
illustrates the effects of self-shielding on the temperature 
structure of the gas more explicitly: the self-shielded gas 
with n > 0.01 cm~'^ is generally cooler (with T < 10'' K) 
when its evolution is consistently modeled. This simply 
results from the suppression of artificial photoheating 
and the enhancement of the cooling function in CIE. At 
these temperatures, the gas radiates very inefficiently 
in Lya, which provides further evidence that the simple 
self-shielding density threshold is not introducing large 
errors: the Lya cooling luminosity versus halo mass 
predicted using prescription 9 (discussed below) is quite 
close to what is obtained by effectively suppressing 
the Lya emission from all the self-shielded gas, as in 
prescriptions 2 and 3 in which the self-shielded gas 
is identified using a ray tracing method and does not 
assume a particular density threshold. We have also run 
a simulation of the Al halo with the ionizing background 
completely turned off, so that all the cooling in this case 
originates from gravitational energy and requires no 
self-shielding correction. The cooling luminosity for this 
simulation equals the one obtained with the simulation 
with on-the-fly self-shielding within ^30%, when the 
star-forming regions are identically excised. We are 
therefore confident that our simple density threshold for 
self-shielding yields relatively accurate results. 

The prescriptions for calculating the Lya cooling 
luminosity from the simulations with on-the-fly self- 
shielding approximation are labeled 6—9 in Table 3 and 
the corresponding results are shown in the right panel 
of Figure 2. For the last three (most consistent) cases, 
the filled symbols show the results for the Al halo in 
our zoom in simulation with 8x better mass resolution. 
Prescription 6 (red x's), in which self-shielding is 
modeled but in which we naively sum over multiphase 
particles, demonstrates how the multiphase particles 
alone can produce artificially high cooling luminosities; 
these should therefore always be excluded, or treated 
separately. As an aside, comparison of prescriptions 1, 
5, and 6 indicates that having the high density regions 
turn into multiphase particles limits the amount of 
artificial photoheating by effectively shielding the very 
dense gas. Three prescriptions (7, 8, and 9) correspond 
to physically plausible cases: one with the on-the-fly 
self-shielding approximation and star formation, but 
with multiphase star-forming particles excluded from 
the Lya luminosity sum (9; blue circles); one also with 
the on-the-fly self-shielding approximation, but with the 
star formation model turned off, summed over all the 
particles (7; magenta squares); and the intermediate 
case of summing only the particles with nn < 0.13 cm~'^ 
at which the gas would have become multiphase if the 



star formation model had been on (cyan diamonds; 8). 
All three are physically realistic in the sense that they 
are uncontaminated by either artificial photoheating 
or by the sub-resolution multiphase model. The only 
difference between the three cases is in how the gas with 
riH > 0.13 cm^'^, the density at which the multiphase 
model becomes active if on, is treated. Since the 
multiphase model was calibrated to match the observed 
Kcnnicutt-Schmidt relation (Springel & Hcrnquist 
2003), this threshold density corresponds approximately 
to the density above which stars should start forming 
(although the exact value depends on some physical 
assumptions and may depend on redshift; e.g., Schaye 
2004; Wolfe & Chen 2006). 

In principle, the second approach (prescription 7) 
might seem more accurate since it captures the entire 
cooling. However, as is apparent in the radiative transfer 
results of §4, the extra cooling luminosity is concentrated 
around the accreting galaxy. It is unclear whether this 
central cooling emission would be observable in reality 
for at least two reasons. First, the density of the 
medium and the immediate proximity of the galaxy 
imply that locally produced Lya photons could be 
efficiently destructed by dust (the escape fraction of 
Lya photons from Lyman break galaxies (LBGs), for 
example, covers the entire range ~ 10""^ — 1; e.g., Kornei 
et al. 2010). In itself, this is not necessarily an issue for 
this study in which we focus on a simplified dust-free 
problem, and would be a well posed problem for follow 
up studies in which dust would be included. Since the 
extra cooling luminosity occurs in ISM gas, it should 
however be accompanied by stellar emission which would 
most likely swamp it locally (see §5), and in that case 
is not really cooling luminosity from the cold streams. 
Second, turning off the multiphase model removes ISM 
pressurization, which can lead to catastrophic collapse 
of the gas rich discs and deepen the potential wells, 
allowing extra energy release. Gravitational interactions 
with dark matter clumps might also artificially transfer 
energy to unstable gas discs. Since this prescription 
assumes that all the baryons are in the gaseous compo- 
nent, while in the central galaxies of actual galaxies in 
massive halos a large fraction would be locked in stars, it 
is likely an upper limit. Prescriptions 8 and 9 are more 
conservative as they exclude all cooling emission occur- 
ring at densities n^ > 0.13 cm^'^. We expect these three 
cases to bracket the true cooling luminosity. A better 
understanding of the energy release at disc interfaces and 
within galaxies will likely be required to make more def- 
inite predictions and should be addressed in future work. 

We conclude by noting that some caution is in or- 
der when interpreting the quantitative details of the 
cooling luminosity predictions for the halos at the low 
end of the mass range in Figure 2, since they contain 
relatively few SPH particles (~ 200 for Mh = 10^° 
Mq) and may not be well converged. Note, however, 
that for the realistic prescriptions 7—9, the slope of the 
numerically predicted La — Mh relation agrees well with 
the analytic expectation based on energy conservation 
(dashed lines in the Figure), in this regime where the 
cold mode dominates and where the scaling should 
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apply. 



4. LYa RADIATIVE TRANSFER 



Having described our hydrodynamical simulations and 
the modeling of Lya photon production, we proceed to 
the Lya radiative transport problem. We use a new 
three-dimensional Lya radiative transfer code, aRT, 
described in more detail in Appendix C. To summarize, 
the fields defining the physical state of the gas from a 
simulation are interpolated onto a Cartesian grid placed 
around a halo of interest, as for the post-processing 
ionizing radiative transfer (§3.2). Monte Carlo Lya pho- 
tons are then seeded throughout the gridded volume, 
with a number proportional to the local Lya emissivity 
(§3). The multiple resonant scatters of each Monte 
Carlo photon are simulated until escape. As the photons 
propagate, the Lya image and corresponding spectrum 
in each pixel on the sky are constructed as seen by an 
observer on Earth, taking into account cosmological 
surface brightness dimming. We do not however model 
the effects of IGM filtering in this work (see §4.2), so 
that the results are more properly interpreted as the 
redshifted emission as it emerges from the galactic halos. 

In typical astrophysical situations, the optical depth at 
the center of the Lya resonance is very large, easily 
To > 10^ or more (eq. CI). As a result, a photon 
propagates only a short distance before being absorbed 
by and exciting another neutral hydrogen atom to its 
2p state. Because of the high Einstein A coefficient for 
the 2p — > Is transition, another Lya photon is quickly 
reemitted {A21 ~ 10~^ s). Since the reemission is in 
general in a different direction than the incident photon, 
the propagation can be effectively pictured as a single 
photon being scattered and undergoing a random walk. 
Even if the scattering is coherent in the frame of the 
scattering atom, the motion of the atom combined with 
the redirection of the photon in general results in a small 
shift in the frequency as viewed by an external observer. 
This results in a random walk in frequency space as well. 
The frequency-space random walk plays a crucial role 
in shaping the emergent line profile as a photon tends 
to escape the medium when it finds itself sufficiently 
far from the line center that the optical depth it sees 
is reduced to ~ 1 (unless the medium is so optically 
thick that the photon spatially diffuses out first). The 
transport of Lya radiation is thus very different than 
ordinary lines and a simplified optically thin treatment 
would lead to fundamental errors both in the theoretical 
predictions and in interpreting observations. 

In this section, we present our basic results on 
Lya cooling emission radiative transfer. As the focus of 
this work is to understand the theoretical and numerical 
uncertainties, and the relevant physical effects, we limit 
ourselves to examining a particular example, our Al 
halo at z = 3. 

4.1. Emission Physics 

In §3, we demonstrated how the predicted Lya cooling 
luminosity depends on the treatment of self-shielding 
and of sub-resolution physics. We illustrate how the 
different possible assumptions manifest themselves in 
other observational properties by performing Lya radia- 



tive transfer calculations for some of the prescriptions 
studied above (Table 3). The results are shown in Fig- 
ures 6 (prescriptions 1, 3, and 4) and 7 (prescriptions 5, 
7, and 8). These fiducial calculations use iVph — 10,000 
Monte Carlo photons. For the radiative transfer plots, 
we subtracted the peculiar motion of the central galaxy 
with respect to the simulation box (40 km s^^ toward 
the observer), so that Av = corresponds to the galaxy 
rest frame. As we had found from our study of a 
sample of halos in §3.2, incorrectly treating self-shielding 
(either by assuming an uniform ionizing background, or 
post-processing ionizing radiative transfer with default 
simulation temperatures) or including gas particles 
contaminated by the effective sub-resolution multiphase 
model can lead to order-of-magnitude overestimates of 
the cooling luminosity. However, it is not only the total 
luminosity of a system that can be incorrectly predicted, 
but also its morphology and line spectrum. 

It is easy to understand how different prescriptions 
for seeding the Lya photons can affect the observed 
morphology and spectrum. In terms of morphology, 
different prescriptions seed the photons not only in 
different amounts, but also in different places. For 
instance, including artificially luminous multiphase gas 
particles produces a disproportionally large Lya lu- 
minosity concentrated in the densest central parts 
of the system, approximating a bright point source 
rather than spatially extended emission from the cold 
streams. Seeding the photons in different places also has 
implications for the emergent spectrum: photons seeded 
deep within dense self-shielded regions have a much 
larger HI column density to traverse before escape, and 
so result in more widely separated double peak profiles. 

Some specific points are worth noting. As indi- 
cated by the surface brightness maps and the circularly 
averaged surface brightness profiles, the resulting objects 
are spatially extended by tens of proper kpc, comparable 
to the virial radius of the halo (see also Dijkstra & 
Loeb 2009) and to many of the observed Lya blobs 
(e.g., Matsuda et al. 2004; Saito et al. 2008). This 
spatial extent results from a combination of some of the 
cooling radiation being emitted in the accretion streams 
far from the galaxy, and of spatial diffusion owing to 
resonant scattering. However, the spatial extent is 
strongly dependent on the surface brightness threshold 
and consequently on the luminosity prescription, so 
that even an intrinsically diffuse source could appear 
relatively compact in observations For instance, for 
our optimistic prescription 7 for the Al halo at z = 3 
shown in the second row of Figure 7, only the central 
few kpc would stick out above the surface brightness 
threshold of « 2 x 10^® erg s^^ cm^^ arcsec^^ of the 
narrowband images of Matsuda et al. (2004). In this 
case, the more diffuse cooling halo would be completely 
missed, but would show up over a larger area in deep 
long slit spectra sensitive to lower surface brightnesses 
(e.g.. Ranch et al. 2008). The predicted line spectra are 
double peaked, a common characteristic of Lya radiative 
transfer reflecting the fact that the photons can escape 
the medium either on the blue side or on the red side 
of the line center, where the opacity is typically too 
extreme (e.g., Neufeld 1990; Zheng & Miralda-Escude 
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Fig. 6. — Dependence of the Lya properties of the Al system at 2 = 3 for prescriptions 1, 3, and 4 for the state of the dense gas. The 
left column shows the observed surface brightness distribution (smoothed with a Gaussian of FWHM = 1"), the middle column shows 
the corresponding circularly averaged surface brightness profile, and the right column shows the line spectrum integrated within the virial 



radius of the halo, indicated by the dashed circles. The surface brightness contours correspond to 10 



10" 



and 10" 



erg s cm 



arcsec"^. The logarithm of the total apparent line luminosity (in erg s~^) within the virial radius is indicated at the top right corner of 
each surface brightness panel. Top: (prescr. 1) Standard hydrodynamical simulation, with uniform ionizing background and a multiphase 
model for star formation, naively integrated over all the particles. Middle: (prescr. 3) Same simulation, but post-processed with ionizing 
radiative transfer (§3.2.1). The self-shielded regions are regions assumed to be in CIE with the temperature set to Tcie = 10,000 K. 
Bottom: (prescr. 4) Same but with self-shielded regions assumed to be in CIE with the temperature set to Tqje = 15, 000 K. 



2002a; Dijkstra et al. 2006; Verhamme et al. 2006). 
In most cases (but not all, reflecting the effects of the 
complex geometry in the different prescriptions), the 
blue peak is slightly more pronounced than the red peak. 
This is a signature of systematic infall in the problem 
at hand, in which the velocity gradients tend to smear 
the line opacity on the red side of the Lya line, making 
it easier for the photons to escape on the blue side. 
Almost all the star-formation powered Lya emitters 
(Shapley et al. 2003; Steidel et al. 2010), as weU as 
many Lya blobs and fainter analogues (Matsuda et al. 
2006; Saito et al. 2008; Ranch et al. 2008), instead 
show dominant red peaks indicative of outflows. We 
do not see this phenomenon here simply because we 
have not modeled galactic winds in our simulations to 
simplify the physical problem. Unlike the thin accretion 



streams (with small covering factor) that produce only 
slightly stronger blue peaks, outflows (with order unity 
covering factor) are expected to boost the red peak 
more drastically (e.g., Dijkstra et al. 2006; Verhamme 
et al. 2006). Intergalactic absorption, neglected here, 
would also tend to preferentially suppress the blue 
peak. In future work, we will include outflows and IGM 
filtering, which should provide a better match to the 
observational data. 

In §3, we had identifled three physically plausible 
prescriptions for calculating the Lya cooling luminosity 
(7, 8, and 9); the radiative transfer results for 7 and 9 
are shown in the second and third rows of Figure 7. The 
surface brightness maps now make it clear that the extra 
luminosity obtained using prescription 7, i.e. when using 
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Fig. 7. — Same as Figure 6, but for prescriptions 5, 7, and 9. Top: (prescr. 5) Simulation with an uniform ionizing background but no 
star formation, integrated over all the particles. Middle: (prescr. 7) Simulation with the on-the-fiy self-shielding approximation and no 
star formation, integrated over all the particles. Bottom: (prescr. 9) Simulation with the on-the-fly self-shielding approximation and star 
formation, but with no Lya luminosity from star-forming particles. Note the different velocity scale in comparison with Figure 6 

the on-the-fly self-shielding approximation but turning 
off the star formation model and summing over the all 
the particles, is concentrated in the densest central parts 
of the system. This also manifests itself in the predicted 
line spectrum, which shows more widely separated peaks 
as a result of the higher column densities of the gas 
through which the bulk of the photons must propagate 
to escape, as discussed in the following section. 

4.2. Radiative Transfer Physics 

Since radiative transfer effects play a crucial role 
in shaping the observational properties, we pause to 
illustrate exactly how important each piece of physics is. 
To do so, we repeat the same fiducial calculations but 
sequentially turn off different physical effects. We only 
repeat the calculations for the plausible prescription 9, 
which suffices to illustrate the physics. The results are 
presented in Figure 8. 



the gas velocities to zero. In this case, the most impor- 
tant change is in the line spectrum, which becomes a 
nearly symmetric double peak. This is easily understood 
in the context of the plane-parallel analytic solution 
derived by Neufeld (1990) (for an adaptation to spherical 
geometry, see Dijkstra et al. 2006) for a monochromatic 
source in an extremely optically thick, static medium. 
Assuming that the source is located at the center of the 
slab and that tq is the line center optical depth from the 
source to the surface, Neufeld (1990) showed that the 
emergent spectrum is a symmetric double peak profile, 
with each peak offset by 



\AVn 



191 km s" 



T 
104 K 



1/6 



\ 1020 en 



1/3 



(5) 



from the center. In dimensionless units in which x — {v — 
i'o)/Ai/£) is the frequency offset in Doppler-broadening 
units, |a:p| = 1.06(aTo)^/'^. For convenience, the corre- 
First, we keep resonant scatters but artificially set sponding offset in observed wavelength units can be ob- 



16 



u 

Q. 
O 



SB (erg s ^ cm 




Faucher-Giguere et al. 
SB Profile (erg s"^ cm^^ "-^) 



Spectrum 



D^ 



le-19 


X 


0.20 


le-20 




0.15 
0.10 
0.05 


le-21 




0.00 



-100 -50 



50 100 




X (proper kpc) 
SB (erg s'^ cm'^ "-^) 



R (proper kpc) 
SB Profile (erg s"^ cm"^ "~^] 



u 

Q. 
O 





41.2 


50 





-50 







^l^ 



-100 -50 



50 100 




X (proper kpc) 
SB (erg s"^ cm"' "-^) 



10 20 30 40 i 

R (proper kpc) 
SB Profile (erg s"^ cm^^ ""^) 



u 

01 
Q. 
O 







41.2 


- 


(3..- 




- 


^~ 


/ 

^ 


f\ 


<\ 





- le-19 l'^) 



0.30 




0.25 
X 


■\ " 


« 0.15 
^ 0.10 


: \ : 


0.05 





-100 -50 



50 100 



X (proper kpc) 



10 20 30 

R (proper kpc) 



40 50 



^O- 



500 500 

Av (km/s) 





Spectrum 


1.0 


/\ rh 


o.a 








0.6 


- / 




\ 


0.4 


- J 


1 1^ 


\ - 


0.2 




\l 





-500 500 

Av (km/s) 





Spectrum 


1.0 


P- 




0.8 


r 


^ 


"S" °'^ 


r 


- 


0.4 


f 


L| 


0.2 


J . 


\ , " 



-500 500 

Av (km/s) 



Fig. 8. — Illustration of the effects of radiative transfer on the morphology and spectrum of the Lya cooling radiation, for prescription 
number 9 (on-the-fly self-shielding approximation with the multiphase star formation model, but with the multiphase particles excluded 
from the luminosity calculation) on the Al system at z = 3. Top: Fiducial calculation, with all radiative transfer effects (repeated from 
Fig. 7 but with different axis scales). Middle: Bulk velocities artificially set zero. The double peaked spectrum is now symmetric since 
the velocity flows no longer break the symmetry between the blue peak and the red peak. The effective line width is set by the physics of 
self-shielding (which determines the HI column densities) and of Lyo radiative transfer. In particular, it has little to do with the global 
properties of the host halo. Bottom: Bulk velocities back on, but resonant scatters artificially turned off. The morphology on the sky is 
now sharper since there is no longer spatial diffusion of the photons. More importantly, the spectrum has lost its double peaked nature 
and its width is set by completely different physics, being in this case representative of the velocity dispersion of the emitting gas. 
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In particular, the emergent line profile and correspond- 
ing effective line width are in this case determined by 
the physics of Lya radiative transfer and of self-shielding 
(which sets the HI column densities) and have little to 
do with the global properties of the host halo. 

Next, we turn velocities back on but artificially turn off 
resonant scattering. This is equivalent to assuming that 
the Lya line were an ordinary, optically thin line. In 
this case, there are two important changes. First, the 
morphology on the sky is slightly sharper since there is 
no longer a spatial diffusion effect arising from photon 



random walks. The spatial diffusion effect appears sub- 
tle here as the cooling emission is intrinsically diffuse. It 
is better illustrated in Figure 9, in which we have kept 
the same physical set up but seeded the Lya photons 
in proportion to the local star formation rate, which 
produces a much more intrinsically compact source. 
To facilitate visual comparison of the morphology with 
the pure cooling cases, we have normalized the star 
formation Lya emission to the cooling luminosity of 
the halo. This halo is however forming stars at a rate 
^ 30 Mq/jt, which could result in up to ~ 3 x 10'*'^ 
erg s~^ of star formation powered Lya emission if the 
escape fraction were unity (Leitherer et al. 1999); in 
this case, it would dominate over the cooling luminosity 
of the halo. Second, the line spectrum is very different 
and has essentially lost its double peaked nature. Most 
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Fig. 9. — Cleaner illustration of the spatial diffusion effect owing to the Lya resonant scatters. Top: Same physical set up as for the 
cooling radiation prescription number 9 (Fig. 8) but with all the Lya photons emitted at the locations of the star-forming particles, in 
number proportional to the star formation rate. To facilitate the visual comparison of the effect on morphology, the total Lya luminosity 
has been normalized to the cooling luminosity for prescription 9. Bottom: Same calculation but with resonant scattering turned off to 
show the intrinsic compactness of the star-forming sources. 



importantly, the spectrum (particularly its width) is 
determined by completely different physics, since it is 
now directly representative of the velocity dispersion 
within the host halo through the Doppler effect. 

An important implication of the radiative transfer 
effects on the observed Lya line width concerns the way 
the masses of extended Lya sources like the Lya blobs 
are estimated. In fact, it is often assumed that the 
Lya line width can be associated with random motion 
within the halo (e.g., Bower et al. 2004; Matsuda et al. 
2006). As could be anticipated from more general 
Lya radiative transfer studies (e.g., Zheng & Miralda- 
Escude 2002a; Dijkstra et al. 2006; Verhamme et al. 
2006) , our results show that neglecting radiative transfer 
effects will tend to overestimate the mass, since they 
alone broaden the Lya line even for completely static 
emitting media. However, the situation is in reality 
more complex since the Lya lines that one actually 
observe on Earth are further filtered by the intervening 
IGM, which can significantly attenuate the Lya line 
(e.g., Zheng et al. 2010a,b; Laursen et al. 2010). At 
z = 3, for example, the diffuse IGM transmits only 
w 67% of photons immediately blueward of Lya and this 
fraction decreases rapidly with increasing redshift (e.g., 
Faucher-Giguere et al. 2008d). Moreover, local matter 
overdensities and systematic infall around massive halos 
(e.g., Faucher-Giguere et al. 2008c) can amplify and 
shift the absorption, whereas galactic winds act to 
"let through" more radiation by redshifting it away 
from the Lya line center (e.g., Santos 2004; Dijkstra 
et al. 2007; Dijkstra & Wyithe 2010). It is therefore in 



general difhcult to relate the observed Lya line alone 
to the properties of the halo producing it, and it is 
not even clear in any given case whether the observed 
line width over— or underestimates the halo velocity 
dispersion. More systematic studies quantifying the 
relation between the observed Lya line properties and 
the host halo mass, extending the type of calculations 
presented in this work, are certainly warranted. We here 
simply caution against taking halo masses estimated 
from the Lya line too literally. This remark applies 
even if the Lya blobs are not predominantly powered 
by cooling radiation since radiative transfer effects will 
necessarily be at play in any Lya source. 

4.3. Convergence 

Most of our radiative transfer calculations employed 
physical fields stored on a Cartesian grid with 256'^ cells 
in a volume of (1 comoving Mpc/h)'^. Since the Lya pho- 
ton trajectories are not tied to the grid points (see Ap- 
pendix C), their increments were however finer than that 
this by a factor of 5. Furthermore, the cell luminosities 
were always calculated directly from the SPH particles 
so that they captured the full gas clumping and were not 
subject to resolution degradation. To verify that our ra- 
diative calculations are robust, we have repeated many 
of them with 512^ grid points in the same volume in- 
stead. As shown in Figure 10 for prescription 9 applied 
to the Al halo at z = 3, the calculations in fact appear 
well converged both in terms of Lya morphology and 
spectrum. 

5. DISCUSSION AND CONCLUSION 
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Fig. 10. — Illustration of the convergence of our radiative transfer results with the number of grid points along each dimension. Left: 
Fiducial cooling luminosity calculation for prescription number 9, with 256^ radiative transfer grid points. Middle: Same, but with 512^ 
radiative transfer grid points in the same volume. Right: Corresponding line spectra within the virial radius. Both the morphology and 
the line spectrum appear well converged with the radiative transfer resolution. 



Motivated by observations of a variety of extended 
Lyof sources at high redshift and by recent studies 
suggesting that cold mode accretion could account for 
the majority of them (e.g., Dijkstra & Loeb 2009), 
we have used hydrodynamical simulations to predict 
the Lya cooling emission of forming galaxies. We 
have introduced a new Lya radiative transfer code, 
aRT, and for the first time applied such a code to the 
particular problem of cooling radiation in cosmological 
hydrodynamical simulations. In this work, which will 
serve as the foundation for follow up studies of the 
astrophysical implications, we have quantified the 
theoretical and numerical uncertainties in predicting 
the properties of Lya cooling radiation. To do so, we 
have considered the simplest physical problem of galaxy 
assembly embedded in the cosmic UV background but 
without feedback processes. We have shown that the 
Lya cooling luminosities, morphologies, and spectra of 
the cold streams are strongly dependent on the treat- 
ment of the ionization and thermal state of the emitting 
gas, and that naive assumptions can produce artificially 
high cooling luminosities. Previous studies have usually 
relied on such assumptions, because most existing hydro- 
dynamical simulations do not self-consistently follow the 
transfer of ionizing radiation (but see Sommer-Larsen 
2006; Laursen & Sommcr-Larscn 2007; Laursen et al. 
2009b), whereas the cold streams in reality self-shield. 
The predicted cooling luminosity is so sensitive to the 
thermal state of the gas principally because of the 
exponential dependence of the coUisional excitation 
coefficient on temperature in the range T '^ lO'' — 10^ 
K characteristic of the cold streams. We have also 
demonstrated that subresolution physics models, when 
not carefully taken into account, can lead to large errors 
as well. 

Having explicitly illustrated the range of uncer- 
tainties, we have made a systematic attempt to converge 
on the correct answer. Our strategy consisted of first 
post-processing our hydrodynamical simulations with 
ionizing radiative transfer to identify the self-shielded 
regions. By inspection, we found that there is a fairly 
well defined total hydrogen density nn ^ 0.01 cm~'^ 
above which the gas self-shields, at least at the redshifts 



z ^ 3 characteristic of the existing observations. We 
have then rerun simulations with exactly the same initial 
conditions and other physical parameters, but with the 
ionizing background turned off in regions exceeding the 
fiducial density threshold nn = 0.01 cm~"^. By turning 
the ionizing background off in dense regions on the fly, 
this dense gas is consistently evolved with modified heat- 
ing and cooling functions, as well as the corresponding 
dynamical response. These simulations provide us with 
our best estimates for the actual Lya cooling luminosity 
of the cold streams, in good agreement with energetic 
considerations. 

For our Lya radiative transfer calculations, we focused 
on a particular halo (Al) of total mass Mh — 2.5 x 10^^ 
M0 at z = 3. We have shown that it is not only the 
integrated Lya luminosity of the system that depends 
strongly on assumptions regarding the thermal state of 
the gas, but also its apparent morphology and spectrum. 
To isolate the role of different physical effects in shaping 
the observational properties, we have sequentially turned 
off separate pieces of physics. Both the Lya resonant 
scatters as well as the bulk velocity structure of the 
system are critical in shaping the emergent spectrum. 
The resultant effective line width in general differs 
strongly from the optically thin expectation owing to 
these radiative transfer effects, and we therefore caution 
against measuring the masses of systems based on the 
Lya line alone. The intrinsically extended nature of 
the cooling emission and the spatial diffusion owing to 
resonant scattering combine to produce objects that are 
spatially extended on the sky, but with an apparent size 
that depends sensitively on the observational surface 
brightness threshold. As a result, only the central few 
kpc of the Al halo Lya cooling emission at z = 3 sticks 
out above the surface brightness threshold ^ 10~^* erg 
s^^ cm~^ arcsec"^ typically achieved to date (Fig. 7). 
If the characteristic size of the Lya emission scales as 

1/3 

R oc Mj^' , then the surface brightness should scale 

as SB ex La/R^ ex La/M^' . For our consistent 
prescriptions 7—9, Figure 2 shows that the high-mass 
luminosity slope varies from about 3/4 to 1 in the most 



optimistic case, i.e. SB ex Mj 



-0.33. 



so that even 



the higher-mass systems are unlikely to show up as 
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sources with extents ^ 100 kpc in existing observations 
from cooling luminosity alone. Fainter sources may 
however well be detectable in deeper observations, such 
as the 100-hr long-slit spectrum reported by Ranch et al. 
(2008). 

Our analysis of a sample of halos from a cosmo- 
logical volume at 2; ~ 3 confirms that (in agreement 
with previous studies; Haiman et al. 2000; Fardal et al. 
2001; Dijkstra & Loeb 2009; Goerdt et al. 2010) coohng 
emission alone can in principle produce luminosities 
La ^ 10**^ — 10^** erg s^^ sufficient to account for 
luminous Lya blobs. This however requires quite 
optimistic assumptions under which most of the energy 
is released close to, or within, the accreting discs, at 
densities sufficient to form stars according to the ob- 
served Kennicutt-Schmidt relation. In those optimistic 
predictions based on prescription 7, star formation was 
artificially turned off (yielding purely gaseous discs), so 
that it is unlikely that this entire cooling emission can be 
realized in reality without being overwhelmed by stellar 
emission.^ Our findings in this respect are at odds with 
the recent simulations of Goerdt et al. (2010), who argue 
that pure cooling radiation can explain the observed 
giant Lya blobs. Our investigation of the sources of 
error in the Lya luminosity predicted from simulations 
(§3) allows us to understand the differences. First, we 
have demonstrated that the predicted Lya luminosity 
is very sensitive to the treatment of self-shielded gas 
and that at z « 3, the characteristic density above 
which gas self-shields from the UVB is nn ~ 0.01 
cm~^. Goerdt et al. (2010) did not explicitly perform 
ionizing radiative transfer and instead assumed a higher 
self-shielding threshold of nn = 0.1 cm^^. This higher 
density threshold translates into an overestimate of the 
Lya luminosity owing to artificial photoheating of gas in 
the density range 0.01 < njj < 0.1 cm~^. Examination 
of the luminosity-weighted T — njj histogram in their 
Figure 7 in fact indicates that the bulk of the Lya lumi- 
nosity they predict originates from this density regime 
(including in clumps that would correspond to satellite 
galaxies, even when the central galaxy is excised), 
whereas Lya emission is strongly suppressed at these 
densities in our approximation to the consistent thermal 
evolution of self-shielded gas, owing to rapid cooling 
below 10"* K (Fig. 5). While there is some uncertainty 
in the precise value of the self-shielding threshold, our 
radiative transfer calculations show that it is clearly 
below 0.1 cm^'^ at the redshifts under consideration (Fig. 
3), and the value nn ~ 0.01 cm~'^ is further supported 
by recent radiative transfer calculations by other groups 
(Nagamine et al. 2010; Aubert & Teyssier 2010). The 
cleanest contrast between our results and those of Go- 
erdt et al. (2010) is provided by their predictions of the 
Lya luminosity originating from the "streams" alone, 
in which they excluded all emission from within 20% 
of the virial radius. Those can be compared, at best, 
with our pessimistic prescriptions 8 and 9, in which only 
gas with rifj > 0.13 cni^"^ (above the subresolution star 

^ In their simulations of LBGs at 2 = 3.6, Laursen et al. 2009a 
find tfiat only ~ 10% of the total Lya emission conies from cooling, 
although wc show in Appendix A that this ratio should depend on 
halo mass and rcdshift. 



formation threshold) is excluded; our predictions are in 
those cases a factor > 10 lower than theirs. It would 
not be fair to compare our most optimistic prescription 
7 with the Goerdt et al. (2010) results with the central 
20% of the virial radius excluded, since most of the 
extra emission under this prescription occurs very close 
to, or inside, the central galaxy. Finally, the simulations 
of Goerdt et al. (2010) included supernova feedback, 
which provides an additional source of energy and could 
also contribute significantly to boosting the cooling 
luminosity, even in the absence of ionizing photons from 
associated local sources (e.g., Taniguchi & Shioya 2000; 
Taniguchi et al. 2001; Ohyama et al. 2003; Mori et al. 
2004). 

It is further notable that many Lya blobs show 
spectral signatures of outflows (dominant red peak), 
at odds with the expectation of a more prominent 
blue peak in the case of pure cooling of infalling gas.^° 
Nonetheless, the release of gravitational potential energy 
through cooling radiation at observationally interesting 
levels is an inevitable prediction of galaxy formation 
in ACDM and the growing body of theoretical work 
supporting the cold mode scenario implies that a sig- 
nificant fraction would come out in Lya. Such sources 
are therefore poised to be routinely detected in future 
surveys, and likely account for at least a subpopulation 
of the more modest extended Lya sources, such as the 
ones detected by Ranch et al. (2008). In fact, at least 
three of the Rauch et al. (2008) sources (#15, 36, and 
37) show spectral signatures suggestive of infall (with 
a more prominent blue peak). Interestingly, GALEX 
non-detections indicate that the bright Lya blobs are 
much rarer at z == 0.8 than at z > 2 (Keel et al. 
2009). This redshift evolution is reminiscent of the 
gradual disappearance of the cold mode at low red- 
shift predicted by theoretical studies (e.g., Birnboim 
& Dekcl 2003; Dekel & Birnboim 2006; Keres et al. 
2005, 2009b). It is important to note that the cold 
mode could play an important role in the existence 
and properties of the Lya blobs even if they are not 
energetically dominated by cooling emission. In fact, the 
presence of cold neutral gas in galactic halos enhances 
the conversion stellar or AGN power into Lya photons, 
and scattering off such gas may be necessary to pro- 
duce the morphological and spectral properties observed. 

Although we have made important progress in ac- 
curately predicting the Lya properties of cold accretion, 
much work remains to be done. For instance, we have 
focused on the simplest physical problem of accretion 
embedded in the cosmic UV background and neglected 
feedback processes. In the actual Universe, galactic 
winds are observed to be ubiquitous at high redshift 
(e.g., Shapley et al. 2003; Steidel et al. 2010) and 
simulations suggest that they are in fact needed to 
reproduce the observed baryonic mass function of 
galaxies (e.g., Keres et al. 2009a; Oppenheimer et al. 
2010). Furthermore, AGN can inject large amounts of 
energy in the surrounding gas (Springel et al. 2005; 
Hopkins & Hernquist 2006; DcBuhr et al. 2009), which 

^^ There are exceptions, for example certain regions of the LAB2 
blob at z = 3.09 (Wilman ct al. 2005). 
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could radiate even after the nucleus has shut off. 
Such processes will modify the kinematic and thermal 
properties of the circumgalactic medium, and therefore 
its Lya signatures. Moreover, they provide additional 
sources of energy that could enhance the total emission. 
The ionizing radiation produced by embedded star 
formation or AGN (e.g., KoUmeier et al. 2010), the 
effects of metals on the cooling function (e.g., Cantalupo 
2010), and the destruction of Lya photons by dust 
(as shown by Laursen et al. 2009b) are also likely to 
be important to reproduce the observed sources, but 
have not been modeled here for simplicity. Since cold 
accretion likely fuels star formation in halo centers, 
at rates similar to the cold gas accretion rates (e.g., 
Keres et al. 2005, 2009b; Faucher-Giguere et al. in 
prep.), gas accretion and star formation are intrinsically 
linked and their roles in Lya emission are difhcult to 
decouple. Follow up studies will build upon the technical 
foundation established in this work and include some of 
these effects. Ultimately, cosmological statistics like the 
luminosity and correlation functions of the Lya sources 
will be predicted. It will also be interesting to con- 
sider complementary observational probes of the cold 
streams, including their Lya polarization (e.g., Rybicki 
& Loeb 1999; Dijkstra & Loeb 2008), their absorption 
signatures, and other emission lines (including from 
metals) . 

On a more basic level, our understanding of the 
physics of the cold mode is incomplete. In particular, 
it is not exactly clear how the cold streams release 
their gravitational energy. At least some simulations 
indicate that they maintain a roughly constant velocity, 
rather than accelerate, as they fall into galactic halos 
(Keres et al. 2005). This implies that the streams would 
continuously radiate the gravitational work done on 
them, perhaps by undergoing a series of weak shocks. 
However, these weak shocks have not been explicitly 
identified in existing simulations. Birnboim & Dekel 
(2003) instead envision a scenario in which the cold 
streams free fall into the halos and release little energy 
before hitting the central disc in a strong shock, which 
may be partially supported by our results which suggest 
that a large fraction of the energy could be released 
near the halo center. Infalling cold gas might also 



release some energy through interactions with halo 
sub-structure and with the lower density, hot halo 
gas. These interactions could involve hydrodynamic 
instabilities below the resolution of current generation 
cosmological hydrodynamic simulations, both SPH and 
AMR. While we have taken the pragmatic point of view 
of predicting the observational signatures implied by 
our current galaxy formation models, it is conceivable 
that the existing simulations are subject to numerical 
limitations. Work is underway to elucidate some of these 
hydrodynamical issues using the new shock-capturing, 
moving mesh code AREPO (Springel 2010). When the 
capability becomes available, it would also be useful to 
perform truly self-consistent ionizing radiative transfer 
on the fly in order to definitively capture the thermal 
evolution of the gas. 
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APPENDIX 

A. ANALYTIC CONSIDERATIONS 

In this section, we outline analytic arguments that motivate the expectation of significant Lya emission from 
galactic gas accretion, at a level comparable to observed Lya blobs. Similar arguments were made by Haiman et al. 
(2000), Dijkstra & Loeb (2009), and Goerdt et al. (2010); their essence is summarized here as a basis to understand 
our numerical results. We also provide an analytic estimate of the amount of power that could be contributed by 
photoionization from the cosmic background. 



We consider gas accreting from the IGM to the bottom of a dark matter halo potential well and ask how 
much energy is available to be radiated. Averaged over the accretion time scale, the gravitational power available 

can be expressed as (i^grav) ~ /cff.Mgas|A<I>gi.av(?'min)|5 where Afgas is the gas mass accretion rate, A<I>gi.av('''min) is 
the potential difference from the IGM to the radius rmin from the bottom of the potential well at which the gas is 
assumed to settle, and /off < 1 is an efficiency factor. The efficiency factor accounts the fraction of the gravitational 
energy that remains in bulk kinetic or thermal form. In addition, only a certain fraction of the energy that is radiated 
comes out in the Lya line. 



Using the fitting formula derived by Neistein & Dekel (2008) for the average halo mass accretion rate, we can 
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write (for the WMAP5 cosmology) 

M,.. « 210 Mo yr-i ( ,f ) ''7^^') '' (^ , (Al) 

where /gas is the fraction of the accreted mass that is gaseous (see also McBride et al. 2009; Fakhouri et al. 2010). For 
a Navarro et al. (1997) halo profile, 

|,^ / M GM20OC In [1 + (ri„inA200c)c] GM20OC c 

|A$grav(rmin)| « — ^ , ,, , , „//-, . ^, ~ — , /-, , „x //-, , . , (A2) 



r„ 



ln(l + c)-c/(l + c) rzooc In (1 + c) - c/(l + c) ' 



iE..„)^-S.S.m"e,,s-^f.„l — Uri ' 5.2 ' ■ <*=" 



where M200C is the halo mass defined as that exceeding 200 times the critical density, r2ooc is the corresponding 
radius, and c is the halo concentration parameter. In the last step, we have assumed cr^i^/r2ooc ^ 1 to simplify the 
expression. In general, M200C ~ M (= MpoF ~ -^isob! see §2.2), but since rigob may exceed r2ooc by as much of 75% 
on cluster scales (White 2002), M200C can be lower than M by a factor of as much as 5 at z ^ 0. For the present crude 
estimate, we will ignore the distinction between M and M200CJ which is much smaller at the high redshifts of interest. 
Combining equations Al and A2, we find 

M _Y' fl + A'-' fU^ \ fc/[ln{l + c)-c/{l + c)] 

This expression is valid at z > 1 (where the cosmological constant can be neglected) and fiducially assumes c — 5, 
approximately the median concentration of all massive halos at z ~ 2 — 3 (e.g., Zhao et al. 2003). 

As an estimate of the sensitivity of this analytic prediction to the shape of the dark matter halos, the calcu- 
lation can be repeated for a Hernquist (1990) profile, 

%^^Ar) = -^, (A4) 

using the relation a — (r2ooc/c)-\/2[ln (1 + c) — c/(l + c)] (e.g., Springel et al. 2005). Wc then obtain the ratio 



(iggrav)"""^"'^* ^ / In(l+C)-C/(1 + C) 
(4rav)^FW V 2 • 

For c = 3, 5, and 10, this ratio takes the values 0.56, 0.69, and 0.86, respectively. Since these differences, at the tens 
of percent level, are subdominant compared to the other sources of uncertainty, we will only plot the prediction for 
the NFW model in this work. 

In the absence of additional sources heating (see below, where we consider the importance of the ionizing 
background), the actual Lya cooling luminosity emergent from a given halo will in general be somewhat lower than 
predicted by equation A3, since some of the gravitational power will be converted to kinetic and thermal forms, rather 
than radiated, and the radiated fraction will not come out entirely in Lya. For halos dominated by the cold mode, 
however, the /off is expected to be significant, perhaps up to ~ 50% if half of the gravitational energy is radiated and 
most of this radiation comes out in Lya. At higher masses, the efficiency will be suppressed as the hot mode becomes 
more important and a higher fraction of the baryons are accreted in the form of stars. 

The analytic model of Dijkstra & Loeb (2009) is more sophisticated than the simple energetic argument above, as it 
accounts for the fact that the fraction of cold gas in halos diminishes with increasing mass. Since only the cold gas 
radiates efficiently in Lya, this yields a shallower slope for the L^ — Mh relation at high masses, which is in fact in 
better agreement with the simulation results (c.f. Fig. 2). Equation A3 is thus really an upper bound that we do not 
expect to be saturated in this regime. 

In addition to gravitational potential energy, there is an additional source of power that can be converted into 
Lya photons even in the case of pure accretion, without feedback from stars or AGN: the cosmic ionizing background. 
An upper bound for how much power the ionizing background can contribute to Lya cooling within a dark matter 
halo can be estimated as the total inward flux of energy from ionizing photons across a sphere of a virial radius. A 
more realistic estimate is obtained by multiplying this quantity by a factor /cov accounting for the fact that only a 
small fraction of these ionizing photons are actually absorbed within the virial shell: 

where i^ion — ''^ J dvJu and J^ is the specific intensity of the ionizing background, assumed homogeneous and 
isotropic, which is a good assumption just above the hydrogen ionization edge at vyj (e.g., Faucher-Giguere et al. 
2009). Since the background spectrum has a strong absorption edge at Av^a owing to intergalactic Hell absorption, 
we can optimistically take J^, « J^.^-^ for v G [vbi, 4i^hi], and J^, — Q beyond Av^a- Then: 

Eion ~ 12Tr^rli^umJiynifcov, (A7) 
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For the hydrogen photoionization rate F ~ 0.6 x 10^^^ s^^ measured from the Lya forest at z = 3 (Faucher-Giguere 
et ah 2008b), this spectrum imphes J^^^ = 1.5 x 10^^^ erg s^^ cm^^ Hz^^. Expressing in terms of halo mass, this 
gives 

"-' " '-' X '»" "'' =- (lO^)"" (.Th)' ( i.5x1O-^^c4T-c„-^H.-0 (si) ■ <*«> 

Comparison of equations A3 and A8 suggests that photoionization from the cosmic background can account for only 
about a percent of the maximum power available from the release of gravitational potential energy at the fiducial halo 
mass Mh = 10^^ M©. Owing to the different mass dependences, and because our numerical results indicate that the 
gravitational power upper bound can be far from saturated in some plausible cases (§3.2), we however cannot exclude 
that the ionizing background contributes a significant fraction of the Lya cooling emission in some regimes. This 
contribution is included in our simulations, though, and while it is difficult to disentangle from gravitational energy, 
circumstantial evidence suggests that it does not dominate. Indeed, in our most realistic simulations with a consistent 
self-shielding approximation the gravitational power upper bound is never systematically violated (see the right panel 
of Fig. 2), and a simulation of our Al halo with the ionizing background completely turned off yields a nearly equal 
cooling luminosity (within 30%) at z = 3 as our prescription 9 when excluding the multiphase star forming regions in 
the same manner. 

Although we intentionally ignore Lya photons produced by star formation throughout most of this work, it is 
interesting to analytically estimate the photoionization stellar Lya emission relative to pure cooling under simplified 
assumptions. Numerical simulations indicate that, in the absence of strong feedback, the star formation rate in 
high-redshift halos closely follows the cold gas accretion rate (Keres et al. 2005, 2009b; Faucher-Giguere et al. in 
prep.), i.e. SFR « Afcoid, except for modest . Stellar Lya emission results from the conversion of the ionizing 
radiation from young stars via absorption and recombination by the surrounding medium, and therefore scales with 
the star formation rate, with a proportionality constant that depends on the initial mass function. For standard 
assumptions, L^^ = lO''^ erg s'^ [SFR/(Mq yr'^)] (e.g., Leitherer et al. 1999). If we assume that Mcoid ~ Mg^s, " 
then we can again make use of the average fitting formula in equation Al to derive the average stellar Lya luminosity 
as a function of halo mass and redshift: 

<Lf ) . 2.1 X 10« e,s .-'/... {J^y {^'f {^^) . (A9) 

where the fa,csc factor quantifies the fraction of Lya photons that avoid destruction by dust. We can then combine 
this with equation A3 to find the characteristic ratio of stellar to cooling Lya: 

(Lf) ffa,csc\f M Y'^^fl + zy^ 5.2 

JS^)""^-^^[1^)[W^M;^) [-^J c/[ln(l + c)-c/(l + c)]- ^ ^^> 

The main point to take away is that it is a priori difficult to predict whether stellar Lya emission should dominate 
over cooling Lya, since the ratio depends on a number of uncertain parameters, including fa.csc a-nd /off. Furthermore, 
the ratio depends on both halo mass and redshift, so that theoretical predictions based on model galaxies at fixed M 
and z (e.g., Laursen et al. 2009a, b) cannot be universally applied. The actual ratio of stellar ionization to cooling 
Lya is likely to be somewhat suppressed relative to the simple scaling in equation AlO because galactic winds can 
remove a large fraction of the gas from the star-forming reservoir (e.g., Murray et al. 2005; Oppenheimer et al. 2010), 
therefore suppressing the Lya emission from star formation ionization, without significantly affecting the accretion of 
cold material and hence the cooling Lya. 

B. IONIZING RADIATIVE TRANSFER 

The opacity in the Lya line (eq. CI) and the Lya emissivity (eq. 1 and 2) depend on the ionization state of 
the hydrogen atoms. The SPH code GADGET calculates the ionization of gas elements assuming photoionization 
equilibrium with a uniform UV background (e.g., Haardt & Madau 1996; Faucher-Giguere et al. 2009) (see §2.3). 
This optically thin treatment misses the effects of self-shielding in dense regions (§2.3). To correct for this, we use 
a post-processing method to solve the radiative transfer problem on the gas distribution provided by GADGET to 
obtain more realistic ionization fractions and to identify the regions that self-shield from external radiation. The basic 
assumption is that the gas dynamics does not significantly depend on the radiative transfer effects missed in the SPH 
calculation. However, re-simulations are also used to approximate the dynamical and thermodynamical response to 
self-shielding (§3.2.2). 

We use a ray tracing algorithm on the same grid as the Lya radiative transfer (Appendix C). Specifically, we 
send a normal ray from each cell on the surface of the radiative transfer volume along the inward direction (x+, x— , 

^^ This assumption is increasingly violated at high halo masses, where hot halo gas becomes more prevalent, but since the cooling 
Lya emission also scales with Meold instead of Mgas to first order, the errors made should approximately cancel when taking the ratio of 
stellar to cooling Lya below. 
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J/+, y~, z+, or z~) and iteratively solve for the equilibrium ionization structure. For this calculation, we take the 
gas temperatures provided by the SPH calculation and assume that all the helium is in the form of Helll. These 
are good approximations for the gas outside of self-shielded regions and they therefore correctly capture the onset 
of self-shielding. In principle, the distribution of helium ionization affects the predicted Lya luminosity through 
the abundance of free electrons. The exact ionization state of helium can however matter significantly only inside 
self-shielded gas, since the free electron number density is dominated by hydrogen elsewhere. Assuming that all the 
helium is in Helll overestimates the number of free electrons and hence the recombination and collisional excitation 
rates. However, the overestimation is only significant for temperatures T < 10^ K, at which the self-shielded gas is 
found to contribute negligibly to the total Lya emission (§3.2). During the re-simulations, turning off the ionizing 
background in dense gas results in most of the helium being in the form of Hel in self-shielded gas, as should be the 
case since Hel traces HI for realistic ionizing spectra. Errors associated with the helium ionization state are therefore 
expected to be small in all cases. The rays are assumed to originate from a diffuse ionizing background with hydrogen 
photoionization rate Fjjj^. For simplicity, we only update the hydrogen ionization fractions and only keep track of 
the photoionization rate FHi.i along each ray, rather than the full ionizing spectrum. This is a fair approximation 
because of the small thickness of the self-shielding layer, in which the gas transitions from almost fully ionized to 
almost completely neutral. Indeed, the thickness of this layer is approximately equal to the mean free path of ionizing 
photons within it: 

A/,, ^ AU^ ^— — « 5 PC ( .."™_3 ) "' , (Bl) 

nm<ym{i'm) VO.Ol cm 3/ 

where cthi is the hydrogen photoionization cross section and i/jji is its ionization frequency, corresponding to 13.6 eV. 
This mean free path is generally much smaller than the spatial resolution of both our SPH simulations (Table 2) and 
our radiative transfer grids. One circumstance in which the self-shielding layer could be significantly thicker would 
be the case in which the illumination is dominated by a nearby quasar (e.g., Cantalupo et al. 2008; Kollmeier et al. 
2010); we do not attempt to model this here, although it would be interesting to explore in the future. 

Along each ray i, the optical depth from the background is calculated assuming that all the ionizing photons 
effectively have a frequency i^ui, Ti = J dlnmam^Viii), and the photoionization rate is correspondingly attenuated 

as Ffji.i = Phi^^^^' • The equilibrium solution is obtained iteratively as follows. Trial values for the ionized fraction 
of hydrogen in each cell are first arbitrarily chosen. The optical depth from the background to each cell, along each 
normal direction, is then computed. The effective photoionization rate within the cell is taken to be the average of 
the attenuated value from each direction. Phi = (1/6) X^i Tni.i- Using this photoionization rate, an updated hydrogen 
fraction is calculated using the balance equation 3. For this step, the electron number density n^ is taken equal to its 
value after the previous update. The procedure is repeated until the ionization fraction has converged in each cell. 
The convergence criterion is set to one part in 100. 

In reality, the gas temperatures provided by the basic hydrodynamical simulations are in error within self- 
shielded regions. We thus store the mean attenuation from the background, (e^^*), in each cell as an indicator of 
self-shielding. Different prescriptions for the self-shielded gas and re-simulation procedures are explored in §3.2. 

C. LYa RADIATIVE TRANSFER: aRT 

In this Appendix, we describe the new three-dimensional Lya radiative transfer code, called ^^aRT" , developed for 
this project. 

Basic Physics and Algorithm 

The basic algorithm is the standard Monte Carlo method, which follows the scatters of a prescribed number, 
A^ph, of Lya photons and statistically estimates the emergent morphology and line profile from them (e.g., Zheng & 
Miralda-Escude 2002a; Cantalupo et al. 2005; Dijkstra et al. 2006; Verhamme et al. 2006; Tasitsiomi 2006a; Laursen 
& Sommer-Larsen 2007; Laursen et al. 2009a). The principal improvements over many of these codes are that it is 
fully three-dimensional and its implementation is parallel to take advantage of distributed computing resources. The 
code therefore scales well to large problems with complex geometries and so can be run on the outputs of large-scale 
hydrodynamical simulations. It produces Lya images, spectra, as well as integral field images (spectra as a function 
of position on the sky). As it was independently developed, its results provide useful checks of existing numerical 
calculations. It has been demonstrated (Tasitsiomi 2006a; Laursen et al. 2009a) that the Monte Carlo algorithm can 
be combined with mesh refinement techniques relatively straightforwardly, which is particularly valuable for problems 
of large dynamic range, for instance when attempting to capture the details of scattering within the clumpy ISM of 
galaxies at the same time as the large scale cosmological radiative transport. We plan to build on the infrastructure 
developed in this work and implement adaptive refinement in future versions of the code. Here, we briefly outline the 
basic physics and methodology of the current implementation. Our notation follows that of Dijkstra et al. (2006). 

There are three relevant reference frames: the frame of the observer, the fluid frame, and the frame of a par- 
ticular atom, which may have thermal motion with respect to the fluid frame. The observer is assumed to be located 
at z — +00 and to be stationary with respect to the simulation volume, apart from a possible Hubble flow that 
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redshifts all the photons from the simulation volume by the same factor. We assume that the simulation volume is 
sufficiently small that its volume is essentially at a single redshift. A quantity Q in the observer's frame is denoted 
Q' in the fluid frame and Q in the frame of the atom. 

The Monte Carlo photons are injected with a probability per unit volume proportional to the local Lya emis- 
sivity (§3). The photons are assumed to have energy exactly at the line center in the frame moving with the fluid and 
the emission is taken to be isotropic. The natural and thermal broadening of the emission line have a negligible impact 
on the results, as the resonant scatters rapidly erase their memory. The Lya photons are then each transported by 
simulating their resonant scatters. 

The Lya optical depth through a HI column density A^hi for a photon at the line center is given by 



ro = A^Hi.Ly.(-o) ^ ^m/.^^^^ « 8.3 x 10^ [ ,^^^-^^_, ) [j^^J^) , (d) 

where crj^yQ, is the Lya cross section, vq is the Lya line center, /12 is the oscillator strength of the transition, e is the 

charge of the electron, rrie is its mass, c is the speed of light, and Aj^d — wth^'o/c, with uth = {'^kBT/nip)^''^ . It is 
convenient to locally define a dimcnsionless frequency x = (i/ — t'o)/^^D- The optical depth for photons of arbitrary 
energy is then 

r. = roH{a, ,) ^ r,- j _^^-—^^^-^, (C2) 

where a = ^2i/(47rAt'D) is the Voigt parameter and A21 is the Lya Einstein A coefficient. 

In the absence of perturbations, the absorption of a Lya photon is quickly followed (after A21 ^ 10^'' s) by 
the reemission of another Lya photon of the same energy. Owing to the motion of the atom, the scatter is however 
not coherent in the observer's frame. Defining Va = Vbuik + vth to be the net velocity of the scattering atom, the 
frequency of a photon of incident frequency Xi after scattering is 

^ — , va.(k.-k.) ,g^^^.^^_,^^o(r-^Y), (C3) 



where k^ and k^ are unit vectors in the directions of the incident and outgoing photons, respectively (e.g., Dijkstra 
et al. 2006). The Hubble expansion is modeled by including a component vh = H{z)r, where r is the displacement 
vector from the center of the system and H{z) is the Hubble parameter, in the bulk velocity term. The g "recoil" term 
accounts for the average transfer of momentum from the photon to the atom during the scatter and can be written as 
g = h/\v-o/2k-QT (Field 1959). The effect of recoil is usually small (Adams 1971) and is ignored in our calculations. 

Because the optical depth near the Lya line center is very large in astrophysical conditions, a Lya photon 
typically travels only a short distance before being scattered. The numerous scatters result in a random walk in space. 
At the same time, the scatters cause the photon to oscillate in frequency and it usually escapes the medium during 
an excursion far into a damping wing of the line, where the optical depth is greatly reduced (e.g., Zanstra 1949; 
Unno 1952; Field 1959; Osterbrock 1962; Adams 1972; Harrington 1973; Neufeld 1990; Gould & Weinberg 1996). As 
a result, the emergent Lya line proflle is heavily affected by the medium through which it propagates and therefore 
provides a probe of the latter. 

To simulate the scatters of the Monte Carlo photons, we at each step randomly pick the optical depth r be- 
fore the next scatter. By deflnition, r has a PDF P{t) oc e""^. Given the frequency and propagation direction of the 
photon, we integrate through the medium until an optical depth t has been reached; this defines the position of the 
next scatter. The frequency seen by a fluid element along the way in general differs from that seen by the observer 
and is given by the Doppler shift formula x[ = Xi — ki • Vbuik/t'th- The frequency of the outgoing photon is then 
calculated using equation C3, which requires picking ko and vth- 

The direction of the outgoing photon is picked from a dipolar phase function, P(ki • kg) ex 1 + (k^ • kg)^. 
Owing to quantum mechanical effects, this Rayleigh phase function is strictly valid only for scatters with Xi > 0.2 
(see Dijkstra & Loeb 2008, and references therein), but the multiple scatters act to quickly randomize propagation 
directions so that the precise phase function is unimportant. An isotropic phase function, P(ki • ko) oc const, gives 
identical results in practically all astrophysical conditions, aside for polarization, which we neglect here. 

The first part of the scattering atom velocity is simply the local bulk velocity of the fluid, Vbuik- To this, a 
thermal velocity Vth must be added. In a frame moving with the fluid, but oriented such that the z' axis is parallel 
to ki, the thermal velocity can be expressed as wth(uj^i, Wj^2' ^11)' where u'u is the normalized component parallel to 
the direction of propagation, and u'^-^ ^ '•^^^ perpendicular components. The probability of scattering off an atom with 
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parallel component u',, is the product of the one-dimensional thermal distribution and the absorption cross section, 

while u'j_i 2 are simply thermally distributed since the probability of absorption is independent of these normal 
components. 

The procedure is repeated for each scatter of each Monte Carlo until it escapes the computational volume. In 
this work, we assume that all photons escape the medium but the algorithm can be straightforwardly extended to 
model photon destruction (e.g., by dust) by specifying a destruction probability at each integration step or scattering 
event. 

Numerical Aspects 

Having described the essential aspects of the algorithm, we now elaborate on numerical aspects of its implementation 
relevant to its accuracy and performance. 

The generation of random numbers is a key element of the Monte Carlo method. In most instances, these are 
straightforwardly generated using either a transformation method or a rejection method (e.g.. Press et al. 1992). An 
important exception is the generation of u'.., whose PDF (eq. C4) is not analytically integrable and for which there 
is not an obvious efficient bounding function to use in the rejection method. Since a variate from this distribution is 
generated at each scatter, it is important to have a fast implementation. We generalize the rejection method described 
in the appendix of Zheng & Miralda-Escude (2002b) to use three critical velocities (in their notation, mi and U2 in 
addition to uq) in order to minimize the number of rejections. 

Although the distance between neighboring grid points is Al = L/Np, photon trajectories are traced exactly 
in the sense that the photon can be located at any coordinate {x, y, z) in the box along its ray. This can greatly 
improve the accuracy of the calculations possible with limited computational resources. For instance, a large and/or 
very optically thick medium may be well described by a relatively coarse grid if it is fairly homogeneous. However, 
because of the short mean free path of the Lya photons, the radiative transfer would be poorly captured if the 
position of the photon were restricted to widely separated grid points. The integration along rays is still done in finite 
steps of length AZjnt, but this parameter can be specified arbitrarily and in particular can be ^ AL Since the photon 
position is in general between the grid points on which physical properties are stored, these must be interpolated 
onto the ray. Our code can use either nearest grid point (NGP) or cloud- in-cell (CIC) interpolation (e.g., Hockney 
& Eastwood 1988). While CIC is more accurate, NGP produces similar results but requires a factor of several fewer 
operations per integration step. We use NGP interpolation in the calculations presented in this work. 

Finally, we use an "accelerated scheme" (e.g., Ahn et al. 2002) that greatly speeds up radiative transfer calcu- 
lations by skipping scatters in the core of the Lya line, during which the photon is nearly stationary in space. 
Specifically, we define a critical frequency Xcrit and say that a photon with \x'\ < Xcrit is in the core of the line. 
Photons in the line core are then forced into the wings by allowing them to scatter only off high-velocity atoms. 
This is achieved by truncating the thermal Gaussian distribution from which u±i and u±2 are drawn to zero for 
\u\ < a:;crit- Skipping core scatters effectively sets the mean free path of core photons to zero. In the limit Xcrit -^ 0, 
the algorithm becomes exact and we can vary Xcrit to verify the convergence of our calculations. In practice, we find 
that Xcrit = 3 generally speeds up the calculations (by as much as orders of magnitude in the very optically thick 
regime) while producing faithful results. We use this value as default in this work. Following Tasitsiomi (2006b, a) 
and Laursen et al. (2009a), we further speed up the code by making use of the fact that in extremely optically thick 
radiative transfer cells (which, in the discretizcd numerical scheme, are internally static and uniform), the solution to 
the radiative transfer problem is well approximated by a rescaled version of the analytic solution obtained by Ncufcld 
(1990). Specifically, the explicit scatters in cells with aro ^ 10'^ can be skipped by choosing the emerging frequency 
according the known analytic solution. 

Visualization 

The principal application of our code is to produce images and spectra that can be compared to observations. To do 
so, the code keeps track of a three-dimensional "integral field" array, of dimensions N^ x N^. The first two dimensions 
correspond to the projected position on the sky as viewed by the observer, or pixels. The third dimension divides 
each pixel into Ni, frequency bins between specified boundaries. From the full integral field array, with a spectrum at 
each pixel, a Lya image can be then be constructed by summing over frequency bins within each pixel. A spectrum 
toward any particular direction, or along an arbitrarily placed slit, can be produced by considering only the relevant 
pixels. 

To construct the integral field array, a number proportional to the probability that the outgoing photon es- 
capes directly toward the observer is added to the corresponding pixel and frequency bin at each scatter. Since the 
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Fig. CI. — Redistribution function test: PDF of the outgoing dimensionless frequency x as a function of the incident xq, for xq = 
0, 1, 2, 3, 4, and 5 from left to right in each panel. The histograms show the PDFs estimated by our Monte Carlo code and the smooth 
curves show the corresponding analytic solutions. The left panel shows the case of an isotropic redistribution function and the right panel 
shows the case of a dipole redistribution function. 

observer is assumed to be located at z = oo, we take this numerical increment to be 



IGtt 



[1 + (k, • zf]e-'^ 



(C5) 



The factor in the square brackets accounts for the probability that the photon is reemitted toward the observer and 
e~'^'' is the probability of escape in a single flight, where Tz is the optical depth to z = oo, given the frequency that the 
photon would have had if it had been scattered in that direction. This procedure (see also Zheng & Miralda-Escude 
2002a; Cantalupo et al. 2005) is necessary because only an infinitesimal fraction of the Monte Carlo photons actually 
emerge exactly in the z — oo direction and is termed the "next-event estimator" (Dupree & Fraley 2002). To convert 
the dimensionless counter to observed physical intensity, we multiply by LQ/[A^ph(A^)^(l + z^], where La is the total 
Lya luminosity and the z factor accounts for cosmological surface brightness dimming. We have verified that this 
next-event estimator conserves the number of photons in isotropic cases, i.e. that the luminosity inferred from the 
resultant image is equal to the luminosity of the central source, from which only a small fraction of the emitted photons 
directly reach the observer. 
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Fig. C2. — Static sphere test; Emergent line profile for a monochromatic source at the center of a uniform static sphere for optical 
depths at the line center tq = 10^, 10®, and 10*. The histograms show the numerical solutions obtained with our code and the smooth 
curves show the corresponding analytical solutions following Dijkstra et al. (2006) (eq. C6). The numerical solution was computed for a 
temperature T = 10 K. 



Test Problems 

In order to be confident in the results obtained from our radiative transfer code, it is important to test it on problems 
with known solutions. Figures CI, C2, and C3 show the results of tests of the aRT code for the following problems: 



Redistribution function. The redistribution function, R{x'g, 



is the PDF of the outgoing photon dimen- 



sionless frequency, x'^, as a function of the incoming photon frequency, x^, in the fluid frame. Although not 
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Fig. C3. — Collapsing and outflowing sphere tests. The left panel shows the numerical results for a collapsing uniform sphere of radius 
R, neutral hydrogen column density A^hi = 2 x 10^" cm~^ from the center to the surface, and temperature T = 20,000 K. In each case, 
the sphere is collapsing with a velocity v = (— r/R)t)niaxf , with Vmax = 20 (dotted), 200 (dashed), and 2,000 km s^^ (solid). The right 
panel shows the same results for an expanding sphere with the same parameters, but reversed velocity field. The results can be compared 
with Figure 7 of Verhamme et al. (2006) for the expanding case. The collapsing case is symmetric under the replacement x —> —x. 

a solution to the radiative transfer problem per se, the Monte Carlo algorithm effectively generates random 
variates from this distribution at each scatter when it picks vth and ko, and applies equation C3. A useful test 
of this key portion of the algorithm is thus to compare the Monte Carlo-generated redistribution function with 
analytically-derived expressions (Unno 1952; Hummer 1962; Lee 1974). Figure CI compares the results of our 
code with the analytic solutions presented by Lee (1974) (for a naturally-broadened intrinsic line) for the cases 
of isotropic and dipole phase functions. To reproduce the redistribution function for core scatters (small x[), the 
accelerated scheme must be turned off (xcrit = 0). 

• Static sphere. An analytic solution to the Lya radiative transfer problem exists for the case of a monochromatic 
point source at the center of an extremely optically thick, uniform, and static sphere of line-center optical depth 
To from the source to the surface: 



J{x) 



r3/2 



\/6aTo 1 1 + cosh [v/27r3/27(|x|3/aT-o)] 



(C6) 



(Dijkstra et al. 2006). This solution is a generalization of the plane-parallel results obtained by Harrington (1973) 
and Neufeld (1990). It differs by a factor of 2tt from the expression in Dijkstra et al. (2006) because it is normalized 
to integrate to unity here. Figure C2 compares the results obtained with our code for tq — 10^, 10®, and 10^ 
and a gas temperature T = 10 K with this solution. 

• Collapsing/outflowing sphere. Velocity fields are crucial to the Lya radiative transfer physics, as Doppler 
shifts modify the perceived optical depth and can open "paths of least resistance" through which photons can 
escape the medium. It is therefore important to test the code in dynamic situations. Unfortunately, few analytic 
results with velocity fields are available. We therefore instead compare our numerical solutions against those 
obtained with other codes. Specifically, we consider the following cases computed by Verhamme et al. (2006) 
(see also Zheng & Miralda-Escude 2002a; Dijkstra et al. 2006): 

- A collapsing uniform sphere of radius R with monochromatic central source: iVui = 2 x 10^" cm^^ 
from the center to the surface, T = 20, 000 K, v = -(r/i?)umaxf with w,„ax = 20, 200, and 2, 000 km s"^ 

- Same, but expanding with v = {r/R)v^^:^r for Wmax = 20, 200, and 2,000 km s^^. 

The numerical solutions obtained with our code shown in Figure C3 agree well with Figure 7 of Verhamme et al. (2006) 
for the expanding case. The collapsing case is symmetric under the replacement x — > —x. 
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